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ABSTRACT 


A  dynamic  model  of  a  multicomponent-multistage  distil- 
lation column  with  variable  holdup  and  boundary  lags  is  developed, 
and  the  programming  for  digital  simulation  is  presented.   The  open 
loop  transient  responses  for  the  mcdel  undergoing  feed  composition 
disturbances  and  reflux  change*  are  studied.   These  responses  are 
then  used  for  synthesizing  transfer  functions  and  closing  the  loop 
for  digital  control  of  the  overhead  product  using  a  variable  reflux 
rate. 
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1.   Introduction. 

Consideration  is  given  here  to  the  dynamic  performance  of  distil- 
lation columns  with  the  ultimate  object  being  column  control  studies 
by  digital  techniques.   This  is  accomplished  by  mathematically  model- 
ing a  multicomponent-multistage  column  and  using  the  model  to  generate 
column  response.   The  column  response  is  then  used  as  an  error  source 
so  that  a  control  algorithm  can  generate  a  corrective  action  to  move 
the  system  back  to  its  desired  control  point.   The  model  performance, 
error  generation  and  control  action  are  all  done  in  a  digital  com- 
puter so  that  in  essence  direct  digital  control  is  being  exercised. 

To  have  realistic  control  action,  model  responses  should  dupli- 
cate actual  plant  responses.   In  order  to  do  this,  it  is  necessary 
to  include  at  the  very  least  column  hydraulics  and  holdup  and i delay 
effects  in  the  column  end  conditions.   In  this  way,  variable  stage 
holdup  is  included  in  the  model  and  the  model  will  respond  to  changes 
in  reflux  flow  which  is  one  of  the  possible  manipulable  variables  in 
this  work.   The  work  of  Peiser  and  Grover  [1]  puts  forth  a  model 
similar  to  the  one  in  this  work. 

It  is  possible  to  generate  model  response  without  these  secondary 
effects.   Models  of  this  type  involve  only  mass  transfer  effects  but 
can  be  quite  complex  if  stage  mixing  and  stage  efficiency  are  incorpo- 
rated.  In  addition  it  is  possible  to  have  interaction  between  hydro- 
dynamics and  mixing  and  efficiency  effects.   In  the  work  presented  here 
it  is  assumed  that  the  stages  are  perfectly  mixed  and  that  the  ef- 
ficiency is  1007o.   This  in  no  way  limits  the  model  developed  since 
allowance  for  both  these  effects  could  easily  be  incorporated  into  the 
model . 

Since  a  model  will  be  of  primary  importance  for  response  gener- 
ation, it  is  important  to  consider  all  possible  physical  processes 
which  can  effect  a  system  response.   Unfortunately,  a  model  to  do 
this  becomes  so  mathematically  unwieldy  as  to  be  useless.   It  has 
become  the  custom  to  make  simplifying  assumptions  which  produce  a 
model  that  is  relatively  easy  to  work  with.   The  sacrifice  made  is 
loss  of  realism  in  responses  produced,  imC  this  is  not  too  important 
for  purposes  of  study  of  general  system  behavior.   However,  to  study 
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system  control  using  models  for  response  generation,  it  becomes 
necessary  to  put  back  in  a  model  some  of  the  "realism"  lost  by  making 
the  "usual  simplifying  assumptions."  The  problem  of  what  to  include 
in  a  model  becomes  a  difficult  one. 

Rosenbrock  [2]  in  his  excellent  paper  on  distillation  column 
control  points  out  that  the  difference  between  observed  and  generated 
system  response  at  small  values  of  elapsed  time  is  mainly  due  to 
secondary  effects  often  not  included  in  dynamic  models.   These  second- 
ary effects  are  the  hydrodynamic  delays  within  a  column  and  the  delay 
and  holdup  effects  in  condenser  and  reboiler  operation  and  external 
mass  transportation.   Unfortunately,  it  is  at  low  values  of  elapsed 
time  that  errors  for  control  action  may  be  generated  and  so  a  model 
to  be  used  for  control  studies  should  include  these  secondary  effects. 

Since  distillation  has  always  been  a  much  used  operation  in 
chemical  engineering,  the  literature  is  voluminous  regarding  steady 
state  analysis  of  multicomponent  systems  and  dynamic  analysis  of 
binary  systems.   The  literature  on  dynamic  analysis  of  multicomponent 
systems  is  somewhat  less  voluminous.   In  particular,  work  on  multi- 
component  systems  wherein  secondary  effects  are  included  is  quite 
limited.   Some  recent  publications  in  this  area  are  the  work  of  Sar- 
gent [3],  that  of  Peiser  and  Grover  [1],  that  of  Moczek  et  al.  [4]  and 
finally  that  of  Anisimov.  [5]  There  are  many  other  papers  (see  for 
instance  [6],  [7],  [8])  besides  these,  but  these  are  representative. 

As  previously  mentioned,  then,  the  work  presented  here  is  the 
development  of  a  multistage-multicomponent  model  which  allows  for 
variable  stage  holdup,,  takes  account  of  liquid  hydraulics  and, 
finally,  allows  for  delays  and  holdup  in  the  condenser  and  reboiler. 
This  model  is  used  to  generate  responses  to  disturbances  in  the  feed 
composition  and  the  reflux  flow.   The  responses  are  used  in  two  ways 
as  follows  (1)  to  abstract  £    -   -.  fer  function  at  any  stage  or  the 
end  of  the  column  and  (2)  to  generate  the  error  needed  for  control 
action  by  means  of  a  numerical  algorithm.   The  transfer  functions 
may  be  used  for  analog  computer  studies  cf  column  control. 

One  of  the  objectives  of  this  work  was  to  see  if  a  numerical 
approach  to  column  control  studies  cc  in  reasonable  time 

using  digital  computers.  As  might  be  expected,  the  answer  to  this 
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is  a  function  of  the  size  of  the  system  studied  and  the  type  of 
control  involved)  but  indications  are  that  work  of  this  type  can  be 
done  in  reasonable  time  at  reasonable  cost.   Another  objective  was 
to  produce  as  accurate  a  model  as  possible  with  the  criterion  being 
duplication  of  actual  plant  response.   Unfortunately,  actual  plant 
data  was  not  available  so  that  model  checking  could  be  done.   Never- 
thelessj  the  model  as  presented  and  used  is  accurate  in  the  sense 
that  at  all  times  during  a  transient,  the  mole  fractions  sum  to  unity 
for  relatively  large  values  of  the  independent  variable  increment. 
This  latter  point  is  considered  in  some  detail  in  the  work  of  Mah 
et  al.  [9]  and  Sargent. [3] 
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2.   Development  of  Model. 


w  ** — 

General  Model 

i  ?      j"1 

F     J 

GENERAL 

STAGE 

(j) 

«— W 

rv    * 

Fu  * 

- — 1-«( 

t 

JU. 

4 

The  above  j  —  stage  (j  is  the  stage  position  index  increasing 
vertically)  shows  all  streams  crossing  the  stage  boundaries.   All 
these  streams  with  the  exception  of  Q  carry  mass  and  energy.   Q,  the 
heat  exchanged  with  the  surroundings,  is  an  energy  stream  only. 

A  general  dynamic  mass  balance  for  component  "i"  at  stage  "j*1  is 


75  (WlVj  +  &  <Vi>l  =  (anpUt  "   ^"^rnass 


(1) 


1  SSF  Wq*    SSF  1.-J 


Eoutput  =  [(Vyp.  +  (LXJ  .]  +  [vsspy      +  LX      ] 


VSSP 


A  general  dynamic  energy  balance  for  the  same  stage  is 


4z   (WTh).  +  -JZ   (W„H)  .  =  (Sinput  -  lout  put) 
d9vLj   d9vVj    x   r         v      ' energy 


(2) 
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Elnput  =   [cH>iml  +   (Lh)j+J  ♦  [fvHfv  +  F^J  +   [vsspHvssF   +  LssphLSSp] 
Output   -    [(VH)  .   +   (Lh)j]  +  Q  +    [Vssp^sp  +  I-SSp\SSp] 


There  will  be  (M»N  +  2)  equations  of  type  (1)  and  (N  +  2)  equa- 
tions of  type  (2).   The  +2   accounts  for  the  boundary  condition  equa- 
tions.  In  addition,  there  will  be  (M»N  +  2)  vapor-liquid  equilibri- 
um relationships  of  the  form 

Kij  "  yi Aj  (3) 

and    (N  +  2)  mole   fraction   constraint  equations 

Plj  =  l  fX±J-l  (4) 

The  equations  as  stated  imply  that  (a)  stage  holdups,  W   and 
W  .,  are  perfectly  mixed,  (b)  equilibrium  between  vapor  and  liquid  is 
instantaneously  realized  at  every  part  of  a  perfectly  mixed  stage, 
(c)  holdup  of  vapor  and  liquid  between  stages  is  negligible,  i.e. 
downcomer  effects  are  not  included,  and  (d)  time  to  attain  fluid 
equilibrium  is  small  compared  with  that  for  mass  transfer. 

Specific  Model 
The  general  model,  while  reasonably  complete,  should  be  modified 
to  reflect  only  the  amount  of  "realism"  required  to  accomplish  the 
purpose  of  the  model.   Otherwise,  the  equations  can  become  unwieldy 
and  needlessly  complex.   With  this  in  mind,  the  following  additional 
assumptions  are  made  about  column  operation: 

(a)  the  column  operates  adiabatically.   (Q  =  0) 

(b)  vapor  holdup  is  negligible  compared  to  liquid  holdup  up  to 
moderate  pressures  and  can  be  neglected. 

(c)  interaction  between  pressure  drop  changes,  liquid  hydraulics 
and  mass  transfer  efficiencies  can  be  neglected. 

These  assumptions  can  be  relaxed  and  the  model  modified  to  in- 
corporate any  of  these  effects  should  it  be  desirable  to  do  so. 

Equations  (1)  and  (2)  are  next  rendered  dimensionless  by  dividing 
through  by  the  feed  rate,  F.   A  dimensionless  time  unit,  t,  is 
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introduced,  defined  by  the  equation 

t  =  (F/W  )  9  (5) 

where  the  time,  measured  in  units  of  tau,  is  equivalent  to  the  number 

of  tray  turnovers  based  on  the  column  feed  rate.   See  Appendix  VI  for 

t,  9  conversion  factor.   Further,  since  it  is  possible  for  W  ,  to  be 

Lj 
different  at  every  stage  and  vary  with  time,  a  particular  W  .  must  be 

LJ 
chosen  for  use  with  the  definition  equation  (5).   A  reference  stage 

is  arbitrarily  chosen  and  the  holdup  value  on  that  stage  at  some  point 
in  time  is  designated  (W  )   and  used  in  equation  (5).   Consistent 
with  the  above  manipulations,  it  is  also  convenient  to  define  a  hold- 
up ratio 

"wj  3  V(Vss  <6) 

Thus,  the  model  equations  (1)  and  (2)  now  become 


dxii         i  r  i 

V   "    Xij   ■  5^  L(EinPUtS    "   Z°UtpUtS)mass   "  VlJ  J 
d^  "  V  ET  [(EinpUt   -   E°UtpUt)energy   "  Vj] 


(7) 


(8) 


Variation  in  liquid  holdup  is  allowed  for  by  using  expressions 
describing  stage  hydraulic  behavior.   Several  detailed  and  very  good 
correlations  describing  such  behavior  have  been  presented  recently, 
i.e.  see  papers  by  Waggoner  [16],  Tetlow  [17]  and  Peiser  and  Grover. 
[1]  A  modified  version  of  the  Peiser  and  Grover  equation  was  adapted 
for  use  here  as  it  provided  an  adequate  but  simple  expression  of  the 
form  W  .  =  :p(L.).   See  Appendix  I  for  derivation.   The  equations  used 


J 


"UJ  =  "PLj  +  atf   \' 


<t^w       »■< 
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<"Lr>SS  [°^  li  J 


S =  *plj  +  fpi.i>iLj  n-p+^itj  frLjJ]  (9a) 

An  overall  mass  balance  yields  the  following  additional  equation 

R,,.  =  (£input  -  Eoutput)  ..  „rtX 

wj       r        r  'mass  overall  (10) 

Several  additional  relations  are  needed  to  account  for  the  coef- 
ficients of  equations  (7)  and  (8)  varying  with  time  if  any  numeri- 
cal solution  of  the  model  equations  is  to  be  attempted.   These  are 
stated  below. 

The  equilibrium  relations  are  assumed  to  be  of  the  form 


lj 


exp"T  h  Bi  +  ciTi 


(ii) 


Normally  K. .  is  a  function  of  temperature,  pressure  and  composition. 
In  equation  (11)  it  is  assumed  that  there  is  no  composition  de- 
pendence and  pressure  drop  through  the  column  is  negligible. 
The  individual  component  enthalpies  are  given  by 


\  -  Vj  +  bi  <12> 


where  the  constants  a,  and  b  are  specific  for  each  component.   This 
same  form  of  temperature  dependence  is  also  assumed  to  hold  for 
individual  component  vapor  enthalpies,  H.,  employing  different 
constants.   The  values  for  these  constants  were  obtained  by  cross- 
plotting  data  from  Maxwell.  [20]  Consistent  with  the  preceding 
relationships  for  enthalpies,  the  enthalpy  of  the  liquid  on  stage 
j,  designated  h.,  is  considered  a  function  of  composition  and 
temperature,  viz. 


hj  =  f(Xij'  V  (13) 


For  other  than  steady  state  conditions,  X,  .  and  T.  are  functions  of 

ij      J 
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time  so  that  the  derivative  of  h  can  be  written 


dh 

dT 


(ah.   dX..]   dh.   dT. 


Assuming  that  h.  can  be  expressed  as  a  linear  combination  of 
individual  component  enthalpies,  namely 

b   h.  =  £h  X  .  (15) 


it  is  possible  to  evaluate  the  partial  derivatives  in  equation  (14) 
and  substitute  back  giving 

V  (?<aiTj +  V  ^J  +[pixij  £1  (16) 

Using  equation  (3)  in  the  form 


?(KiXi)j  =  l  (17) 


and  taking  the  time  derivative  gives 

?(vu  +  xu  51]  =  °  (18> 


where 


dK 

~d9 


f      A.      /dT.v  /dT.v) 


Substituting  (19)  and  (18)  into  (16)  and  converting  to  dimensionless 
time,  t,  gives 
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(EiAi 


h   =  £(a.T.  +  b.)X4.  +  Sa.X.  .1  '  'J  -J  >         (20) 


The  changes  in  composition  from  the  integration  of  the  mass 

a 

balance  equations  (7)  across  At  give  values  of  X.  .  and  X. .  at 

(t  +  At).   From  the  X  . ,  the  T.  at  (t  +  At)  are  obtained  by  a  bubble 

point  calculation.   Thus,  equation  (20)  provides  an  independent 

estimate  of  h  since  all  of  the  terms  are  known  at  (t  +  At). 

J 
Equation  (20)  used  in  conjunction  with  equations  (8),  (9)  and  (10) 

suggests  an  iterative  scheme  to  correct  internal  column  flows  $ 

across  time. 

In  summary,  then,  equations  (4),  (7)  through  (12),  (15),  and 
(20)  are  the  model  equations  used  for  a  numerical  solution  of  the 
general  stage  dynamic  response. 

It  might  be  noted  here  that  the  individual  component  mass  balance 
equation* (7)  can  be  considered  to  be  a  finite  difference  approxi- 
mation to  (M«N  +  2)  coupled,  nonlinear,  partial  differential  equations 
in  X. .,  the  independent  variables  being  time  and  distance  along  the 
column.   These  partial  differential  equations  are  of  second  order  in 
distance  along  the  column  and  first  order  in  time  and  also  contain 
some  first  order  terms  in  distance.   This  would  indicate  that  the 
solutions  do  indeed  describe  traveling  waves  of  composition  changes. 
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3.   Numerical  Solution. 

Preliminary  Discussion 
The  data  and  results  of  this  study  were  obtained  using  a  CDC 
1604  digital  computer  to  simulate  the  column.   This  section  will  treat 
the  numerical  integration  of  the  general  stage  equations.   The  actual 
programs  used  for  various  boundary  conditions  are  given  in  Appendix 

II. 

Before  it  is  possible  to  simulate  a  transient  response,  it  is 
necessary  to  know  the  initial   steady  state  conditions  for  the  nu- 
merical integration.   All  initial  and  final  steady  state  conditions 
were  calculated  using  Program  I  given  by  Hanson  et  al .  [11]  Results 
for  the  unperturbed  and  perturbed  steady  state  conditions  are  tabu- 
lated in  Appendix  VII.   Also,  since  the  constants  in  the  hydraulic 
equations  (9)  contain  specific  column  parameters,  some  column  design 
parameters  for  a  typical  column  must  be  specified.   See  Appendix  III 
for  values  used  in  this  work. 

General  Numerical  Scheme 
A  general  approach  for  obtaining  a  numerical  solution  to  a  system 
of  nonlinear  differential  equations  with  time  varying  coefficients  is 
to  set  up  a  trial  and  error  sequence  for  the  integration  with  correct- 
ive iteration  of  the  coefficients  eventually  forcing  the  equations  to 
be  satisfied.   As  was  mentioned  in  the  section  on  the  specific  model, 
equations  (8)  and  (20)  suggest  an  iterative  scheme  for  correcting 
internal  column  flows  within  each  time  increment.   Following  this 
approach,  the  logic  sequence  used  in  programming  a  numerical  solution 
for  the  response  is  as  follows: 

(a)  one  stage  is  designated  as  the  starting  point  for  beginning 
the  integration. 

(b)  equations  (7)  for  that  stage  are  integrated  across  one  At 
using  the  Runge-Kutta-Gill  algorithm  as  presented  by  Gill.  [10]  The 
coefficients  are  held  constant  at  their  last  known  values.   This 
moves  the  X. .  to  the  point  (t  +  At). 

(c)  a  bubble  point  calculation  using  equations  (3),  (4)  and  (11) 

is  performed  to  establish  a  new  T.. 

J 

(d)  liquid  and  vapor  stream  enthalpies  are  calculated  from 
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equations  (12)  and  (15). 

(e)  equations  (8)  and  (20)  are  evaluated. 

(f)  the  test  h.     /onx  -  h       J  <  e  (an  arbitrary  small 
v  '  I  j  eqn(20)     j  eqn(8)| 

number)  is  now  used  for  convergence  within  a  single  At.   If  the  test 

fails,  for  the  first  iteration  only,  R^   is  incremented  by  an  arbi- 

•  WJ  1 

trary  amount,  AR^  ,  and  R  .  is  evaluated  as  AR./At.    L.  is  calcu- 
lated from  equation  (9),  V.  is  calculated  from  equation  (10),  and 
steps  (b)  through  (f)  are  repeated. 

(g)  for  the  second  and  subsequent  iterations,  the  last  two  known 

values  of  R^  are  used  to  predict  by  the  method  of  "reguli  falsi" 

(21)  that  value  of  R   which  will  cause  the  condition  in  (f)  to  be 

2         *"  J 
satisfied.   When  this  condition  is  met,  it  is  possible  to  proceed 

to  the  next  stage  and  so  on  up  and  down  the  column  until  all  stages 

have  been  integrated. 

(h)  increment  At  and  begin  sequence  again  at  (a).   A  flow  diagram 

of  this  procedure  is  given  in  Appendix  IV. 

Boundary  Conditions 
In  general,  some  control  must  be  exercised  or  the  column  will 
not  operate.   For  initial  investigations  it  was  assumed  that  the  feed 
rate  is  flow  controlled  and  the  feed  is  saturated  liquid  at  its 
bubble  point;   also,  top  product,  bottom  product  and  reflux  stream 
flows  were  set  and  held.   Thus,  any  changes  in  internal  column  flows 
are  absorbed  as  holdup  changes  entirely.   In  a  later  section  when 
control  is  discussed,  these  restrictions  are  modified.   But,  for 
the  present,  these  assumptions  are  satisfactory  for  checking  model 
accuracy  and  looking  at  the  open  loop  response. 

This  was  found  to  be  a  good  approximation  after  comparing  results 
using  equation  (9a)  and  using  an  arbitrary  AR^  .   Use  of  equation  (9a) 
requires  a  more  complex  iteration  scheme  whicft-'is  very  sensitive  to 
changes  in  internal  flows  and  fails  to  converge  except  at  very  small 
increments  of  At. 

2 
In  steps  (f)  and  (g)  several  other  similar  schemes  were  tried 

which  involved  equating  equation  (8)  to  equation  (20)  and  solving 

for  L..   R  .  and  k    .  were  then  calculated  directly  using  the  new 

valueJof  L..   All  sdch  attempts  failed  to  produce  convergence  within 

a  single  ^  At. 
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A  total  condenser  was  used  in  all  simulation  runs,  both  with 
and  without  an  accumulator.   Thus,  reflux  liquid  is  at  its  saturation 
temperature  at  all  times.  When  an  accumulator  was  incorporated,  it 
was  assumed  to  be  an  adiabatic  mixer  undergoing  composition,  temper- 
ature and  holdup  changes  which  in  turn  specify  reflux  stream  con- 
ditions. 

The  reboiler  was  treated  as  a  large  capacity  equilibrium  stage 
with  a  heat  exchanger.   No  heat  exchanger  dynamics  were  included, 
but  can  be  incorporated  at  any  time.   The  reboiler  duty  was  treated 
in  two  ways . 

(a)  Reboiler  duty  (Q  )  was  set  and  held  at  the  new  known  steady 

R 

state  value  during  the  entire  transient  response. 

(b)  Q0  was  considered  a  variable  and  allowed  to  float  within  a  At 

R 

increment  toward  a  new  value  which  would  correspond  to  the  new  equi- 
librium conditions  in  the  reboiler  at  that  point  in  time. 

The  first  situation  was  adopted  for  runs  presented  here.   Physi- 
cally, this  means  that  the  controller  set  point  for  energy  input  to 
the  reboiler  is  not  changed. 

Delays 

A  major  use  of  the  model  response  will  be  to  derive  transfer 
functions  for  control  studies  using  the  analog  computer.   In  the 
many  possible  control  schemes,  it  may  well  be  that  errors  are  detected 
at  small  values  of  elapsed  time.   Any  time  delay  in  the  fluid  trans- 
port at  the  column  ends  will  affect  the  derived  transfer  function,  so 
the  model  should  have  transport  lag  effects  incorporated  in  it.   One 
point  where  these  lags  occur  are  in  the  vapor  flow  from  the  top  stage 
through  the  condenser-accumulator  and  reflux  return  line.   Another  is 
in  the  vapor  flow  in  the  reboiler  return  line. 

To  introduce  delays  at  these  two  points,  a  simple  queing  tech- 
nique was  used.   A  plug  flow  situation  was  assumed  in  which  the 
entering  vapor  stream  variables  (V.,  T.,  H.,  y. .)  are  those  leaving 
the  top  plate  or  reboiler  at  a  particular  increment  in  time  and  that 
these  variables  can  be  moved  across  the  delay  interval  without  mixing. 
Thus,  the  variables  leaving  the  delay  interval  correspond  to  con- 
ditions "n"  At's  ago.  Diagrammatically,  this  can  be  represented  as 
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follows: 


(vj  Tj  hj  V 


At 


ATl 


AT; 


7    I 


At 


(V.   T.    H.    y.  .)    . 
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4.   Numerical  Results 

Column  Response  to  a  Feed  Perturbation 
To  check  the  model,  computer  runs  were  made  with  a  perturbed 
feed  and  variable  end  conditions.   See  Fig.  1  for  physical  system. 
The  feed  entering  plate  4  consists  of: 

Hydrocarbon  Original  Perturbed 

Composition  Composition 

C3  .25  .25 

C,  .25  .30 

4 

C5  .25  .20 

C,  .25  .25 


Variable  boundary  condition  employed: 

w 

1 


(a)  model  with  total  condenser,  constant  R„.>  no  accumulator 


and  no  delays. 

(b)  model  as  in  (a)  with  variable  R  .  added. 

(c)  model  as  in  (b)  with  accumulator  added. 

(d)  model  as  in  (c)  with  condenser  and  reboiler  delays  added. 

2 
Moczek,  Otto  and  Williams   [4]  reported  a  substantial  effect  of 

the  computation  time  increment  used  on  the  accuracy  of  the  transient 

response;   this  was  reflected  to  a  large  extent  in  their  derived 

3 
approximate  transfer  functions.  Rosenbrock  also  discussed  this 

point  earlier.   The  main  idea  is  the  importance  of  using  a  compu- 

For  a  more  complete  discussion  of  this  model,  refer  to  a  paper 
by  Duff in,  J.  H. ,  "Improving  the  Accuracy  of  the  Dynamic  Response 
of  a  Multistage-multicomponent  Column  Model,"  AIChE  Symposium  on 
Systems  and  Process  Control,  56th  National  Meeting,  San  Francisco, 
California,  May,  1965. 

2 
Moczek,  J.  S.,  Otto,  R.  E.  and  Williams,  T.  J.,  "Approximation 

Models  for  the  Dynamic  Response  of  Large  Distillation  Columns," 

AIChE  Symposium  on  Process  Control  and  Applied  Mathematics,  Vol.  61, 

1965,  pp.  136-146. 

3 
Rosenbrock,  H.  H. ,  Brit.  Chem.  Eng.,  Vol.  3,  1958, 

pp.  364-367,  432-435,  491-494. 
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4.   Numerical  Results 

Column  Response  to  a  Feed  Perturbation 
To  check  the  model,  computer  runs  were  made  with  a  perturbed 
feed  and  variable  end  conditions.   See  Fig.  1  for  physical  system, 
The  feed  entering  plate  4  consists  of: 

Hydrocarbon  Original  Perturbed 

Composition  Composition 

C3  .25  .25 

C,  .25  .30 

4 

C  .25  .20 

C,  .25  .25 


Variable  boundary  condition  employed: 

(a)  model  with  total  condenser,  constant  R^.,  no  accumulator 
and  no  delays. 

(b)  model  as  in  (a)  with  variable  R  .  added. 

(c)  model  as  in  (b)  with  accumulator  added. 

(d)  model  as  in  (c)  with  condenser  and  reboiler  delays  added. 

2 
Moczek,  Otto  and  Williams   [4]  reported  a  substantial  effect  of 

the  computation  time  increment  used  on  the  accuracy  of  the  transient 

response;   this  was  reflected  to  a  large  extent  in  their  derived 

3 
approximate  transfer  functions.  Rosenbrock  also  discussed  this 

point  earlier.   The  main  idea  is  the  importance  of  using  a  compu- 

For  a  more  complete  discussion  of  this  model,  refer  to  a  paper 
by  Duff in,  J.  H. ,  "Improving  the  Accuracy  of  the  Dynamic  Response 
of  a  Multistage-multicomponent  Column  Model,"  AIChE  Symposium  on 
Systems  and  Process  Control,  56th  National  Meeting,  San  Francisco, 
California,  May,  1965. 

2 
Moczek,  J.  S.,  Otto,  R.  E.  and  Williams,  T.  J.,  "Approximation 

Models  for  the  Dynamic  Response  of  Large  Distillation  Columns," 

AIChE  Symposium  on  Process  Control  and  Applied  Mathematics,  Vol.  61, 

1965,  pp.  136-146. 

3 

Rosenbrock,  H.  H. ,  Brit.  Chem.  Eng.,  Vol.  3,  1958, 

pp.  364-367,  432-435,  491-494. 
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tation  interval  of  a  certain  minimum  size  even  though  computer  time 

is  correspondingly  increased.   The  interval  may  be  varied  over  the 
course  of  a  computation  to  give  desired  accuracy  and  speed  up  solution 
time  at  larger  values  of  elapsed  time. 

For  this  study,  the  computation  interval  (At)  was  varied  from 
.01  to  .3  (comparable  to  the  .02  to  3.0  minute  interval  used  by 
Moczek,  Otto  and  Williams)  and  very  little  effect  on  accuracy  was 
noticed.   This  was  attributed  to  the  small  number  of  plates  used 
in  this  model.   See  Fig.  2.   In  Fig.  2  the  short,  solid,  vertical 
lines  indicate  the  spread  of  points  for  At  varying  between  .03  and 
.30;   the  dotted  lines  join  the  data  points  for  At  =  .1.   The  value 
of  At  =  0.3  was  found  to  be  the  upper  limit  for  a  stable  solution. 
A  At  of  .1  was  selected  for  use  in  all  subsequent  runs  although 
computer  time  to  attain  steady  state  was  on  the  order  of  45  minutes  with 
ample  print  out  data  for  all  runs.   No  attempt  was  made  to  vary  At 
during  a  transient  solution. 

Figures  3-8  are  temperature,  mass  transfer,  liquid  flow  and 
holdup  response  curves  representing  several  stages  in  the  column  for 
varying  column  end  conditions.   These  runs  provided  a  check  on  the 
accuracy  of  the  model  in  the  sense  that  the  model  did  move  toward  a  known 
steady  state  with  the  mole  fraction  sums  always  very  close  to  unity 
and  all  time  derivatives  tended  to  zero  within  the  error  limit 
allowed  by  the  iterative  scheme.   Table  1  lists  representative  values 
of  these  accuracy  criteria  over  a  typical  solution  interval. 

The  effect  of  including  tray  hydraulics,  lags  and  delays  should 
be  a  slowing  of  the  mass  transfer  response  with  the  final  steady  state 
being  the  same  for  all  cases.  Referring  to  Fig.  5,  the  temperature 
response  of  the  top  and  bottom  stages,  and  more  specifically  Fig.  3 
and  4jtop  stage  and  feed  stage  composition  response  for  the  C, 
fraction,  it  is  evident  that  this  is  indeed  the  case.   The  responses 
flatten  out  and  the  approach  to  equilibrium  becomes  slower  as  the 
model  is  altered  to  include  additional  effects.   An  enlargement  of 

Computer  run  time  was  substantially  reduced  by  reprogramming 
in  FORTRAN  63  vice  60. 
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the  initial  portion  of  Fig.  3  which  is  shown  in  Fig.  9  revealed  some 
dead  time  except  for  the  constant  holdup  model.  With  only  eight 
stages,  dead  time  was  very  small  and  noticeable  only  three  to  four 
stages  away  from  the  perturbation.   For  all  intents  and  purposes, 
responses  were  immediate  within  the  column. 

Referring  to  Fig.  6  and  7,  the  liquid  flow  responses  on  the  top 
stage,  feed  stage  and  bottom  stage  were  noticeably  different  from 
the  mass  transfer  responses.   The  addition  of  hydraulics,  an  accumu- 
lator and  then  delays  accelerated  initial  changes  in  the  liquid  flows. 
At  larger  values  of  elapsed  t,  the  response  curves  crossed  so  that 
final  approach  to  equilibrium  was  again  slower  as  more  effects  were 
added  to  the  model.   It  can  be  rationalized  that,  while  the  mass 
transfer  responses  show  expected  behavior,  the  lags  introduced  in  the 
condenser-accumulator-reflux  cycle  and  the  speed  of  propagation  of 
disturbance  waves  in  liquid-vapor  flows  resulting  from  the  inclusion 
of  hydraulic  interactions,  cause  changes  in  equilibrium  conditions 
which  badly  upset  the  vapor-liquid  flows  initially.   The  rapid  shift 
in  internal  flows  should  smooth  out  as  the  disturbance  waves  damp  out. 
Also,  the  iterative  scheme  for  movement  of  column  variables  implicit- 
ly alters  vapor-liquid  flows  to  meet  the  convergence  test  within  each 

At.   This  produces  rapid  initial  changes  in  flows  which  gradually 

»    • 
smooth  out  as  the  time  derivatives  (h  ,  R.,.)  pass  through  their 

J   Wj 
maxima  and  become  small.   The  final  approach  to  equilibrium  should  be 

similar  to  the  mass  transfer  responses. 

Representative  effects  of  hydraulic  equations  on  feed  stage 

holdup  are  illustrated  in  Fig.  8.   Table  2  gives  the  magnitude  of 

holdup  changes  for  all  eight  stages.   It  should  be  noted  that  care 

must  be  exercised  in  interpreting  responses  at  very  low  elapsed  t 

values.   Any  mismatch  in  the  steady  and  unsteady  state  model  will 

cause  initial  disturbances  which  must  be  allowed  to  die  out  before 

any  desired  disturbance;?  is  imposed  on  the  system. 

Column  Control  Studies 
Control  in  the  digital  sense  was  done  by  letting  the  computer 
correct  the  responses  through  a  control  variable.   This  variable  was 
adjusted  via  a  control  algorithm  based  on  a  generated  error.   The 
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FIGURE   6 
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FIGURE   7 
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FIGURE  9 


control  objective  was  to  maintain  constant  C,  composition  in  the 
overhead  product  by  varying  the  reflux  to  compensate  for  disturbances. 

Accordingly,  the  top  product  response  to  reflux  changes  was 
needed  to  determine  the  proper  algebraic  form  of  the  control  equation. 
Also,  these  open  loop  responses  are  necessary  for  synthesizing  transfer 
functions  for  analog  computer  studies.  Fig.  10  shows  the  actual  model 
response  to  reflux  perturbations.  Fig.  11  is  a  plot  of  reflux  rate 
versus  overhead  C,  composition  at  final  steady  state.   Results  of 
these  runs  indicate  that  operation  at  a  reflux  rate  around  1.0  limits 
the  control  action  severely.   In  the  vicinity  of  a  1.0  reflux  rate, 
the  steady  state  composition  goes  through  a  maximum.   Thus,  reflux 
changes  would  be  ineffective  in  compensating  for  decreases  in  overhead 
composition  at  that  point.   As  long  as  the  perturbations  into  the 
system  produce  increases  in  overhead  G,  composition,  then  control  can 
be  maintained  in  this  region.   By  shifting  the  initial  operating  point 
away  from  1.0  reflux  rate,  fairly  effective  control  can  be  realized. 

It  was  interesting  to  note  also  how  fast  the  reflux  disturbances 
were  propagated  through  the  column.   Fig.  12  shows  that  there  is  a 
small  delay  before  any  change  in  liquid  flow  is  detected  in  the  bottom 
stage  and  a  slightly  longer  delay  before  any  composition  change  occurs. 
This  4s  expected.   It  takes  a  finite  amount  of  time  for  disturbance 
waves  to  travel  the  column  length.   If  the  model  is  performing  properly, 
with  only  eight  stages  the  delay  should  be  slight,  but  detectable. 

For  this  control  study  the  model  contained  the  following: 

(a)  variable  holdup 

(b)  delays  in  condenser  and  reboiler  lines 

(c)  large  capacity  accumulator 

(d)  liquid  level  control  in  accumulator  and  reboiler 

The  computer  was  programmed  to  use  a  proportional  plus  integral 
control  equation  of  the  form  r  =  r0  +  K_e  +  K   f^T 
fa here  the  error,  e,  was  generated  by  comparing  the  instantaneous 
value  of  C,  in  the  accumulator  to  the  initial  starting  value.  An 
arbitrary  error  band  of  .001  in  mole  fraction  was  imposed,  i.e. 
control  action  was  taken  if  jx  -  Xq  J  >  .001.   X  denotes  composition 
of  C,  in  the  accumulator.   Fig.  13  shows  the  results  of  varying  K 
and  K  .   The  values  K  =  K  =  1,  gave  an  underdaraped  response  with  a 
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TABLE  1  .  i 

MOLE      FRACTION    SUM    a      DERIVATIVE      VALUES    ON 

FEED     STAGE      FOR     A    TYPICAL     COMPUTER    SOLUTION 

(Ar=O.I) 


MODEL 

WITH 

1 

HYDR/J 

WLICS 

MODEL 

WITH    HYDRAULICS  8 
ACCUMULATOR 

#OF 
T'S 

i'xi 
i 

EON  8 

EON  13 

i 
Rw 

i 

• 
h 

EON  8 

EON  13 

Rw 

1 

.99602 

-189.9 

-  190.3 

-3 

1.2X10 

.99682 

•189.9 

-190.3 

I.2XI0*3 

10 

1.00008 

-92.2 

-91.8 

-4 
9.5X10 

1.0008 

-92.2 

-91.9 

9.7XI0"4 

00 

.99999 

-35.6 

•35.9 

-4 
3.0X10 

1.00003 

-37.5 

-37.1 

4.9XI04 

100 

1.00000 

-24.3 

-23.8 

-4 
5.0X10 

1.00001 

-24.6 

-24.1 

4X)XI0*4 

900 

.99999 

-3,9 

-3.6 

-4 

I.OXIO 

.99999 

-4.7 

•4.2 

1.8  XIO*4 

1000 

.99999 

-.7 

-.7 

-5 
8.0X10 

1.00000 

-.5 

-.6 

-5 

7.0X10 

♦800 

.00S09 

♦.e 

*•' 

-5 
6.0X10 

.99999 

-.4 

-.4 

-5 

e.oxttr 

TABLE 
INITIAL  a    FINAL 


VALUES     OF   HOLDUP     RATIO 


STAGE  # 


Rw  INITIAL 


Rw  FINAL 


AR, 


REBOILER 


10. 0000 


10.2100 


.2100 


.9859 


9622 


.0037 


9946 


.9971 


.0025 


9991 


1.0087 


.0098 


1.0000 


1.0145 


.0146 


.6442 


8673 


.0231 


8641 


.8913 


0272 


8826 


.9050 


0224 
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long  settling  time.   Increasing  the  proportional  gain  to  five, 
resulted  in  an  overdamped  response  which  settled  out  very  quickly. 
However,  a  K  =  K  =  10  produced  an  unstable  response  with  increasing 
amplitude  of  oscillation.   No  attempts  were  made  to  optimize  the 
response  in  any  way.   This  portion  of  the  work  was  merely  to  ascertain 
if  control  could  be  done  in  a  reasonable  time  using  digital  control. 
Results  were  encouraging  in  this  respect  and  the  model  performed 
adequately  in  all  cases. 
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5.   Transfer  Function  Synthesis. 

Technique  Used 

The  purpose  of  this  section  is  to  outline  the  technique  used 
in  synthesizing  an  approximate  transfer  function  for  the  column  and 
present  the  results  of  one  specific  case. 

The  Laplace  transform  of  X(t), 


■1 


|e"StX(t)dt  (21) 

0 

can  be  transformed  to  the  interval  (0,1)  by  means  of  the  substitution 

-t 
r  =  e 


Then 


X(s)  =  |  r   X(-log  r)dr  (22) 

0 

This  integral  can  be  approximated  by  a  Gaussian  quadrature  formula 
of  order  N  appropriate  to  the  interval  (0,1), 

N  s-1 

X(s)  *  S   v  X(-log  r.)r.S  (23) 

i=l 


s  =  1,2, 


where  (r.)  are  the  roots  of  the  shifted  Legendre  polynomialss 

*     x 
PN  00. 

PN*(X)  =  PN(1-2X) 

and  (w  )  are  the  corresponding  Christoffal  weights. 

Tabulated  values  of  the  (r.)  and  the  (w.)  are  given  for  N  up 

to  15  in  Bellman,  et  al.  [15].  XSee  also  Beliman,  et  al.  [12]  for 
general  approach  to  numerical  evaluation  of  Laplace  transforms. 
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Use  of  equation  (13)    in  its  present  form  restricts  the  X(t) 

values  to  a  very  small  interval.   The  upper  limit  of  the  interval  does 
not  increase  significantly  as  the  order  of  N  is  increased,  but  the 
computational  effort  does.   Therefore,  to  expand  the  X(t)  interval  by 
a  factor,  "a,"  one  makes  use  of  the  multiplicative  property  of  the 
transform  and  equation  (Z3)  becomes 

x(sJa)   =  r      w  X(-a  log  r  )r  8"1  (24) 

i=l 


=  1,2, 


If  the  functional  form  of  X(s)  is  known  approximately  and  it  is 
desired  to  improve  the  estimate  of  parameters  contained  in  X(s),  one 
can  use  a  least  squares  curve  fitting  technique  on  the  values  of 
X(s)  from  equation  (24)  and  the  calculated  values  of  the  approximate 
or  assumed  functional  form  of  X(s),  i.e. 


N     f  2  2  1 

sjl<.)   eqn(4)   -  X(s)assimedJ 


where  S  is  to  be  minimized  with  respect  to  the  parameters  in 

X(s) 

assumed 

Curve  fitting  in  the  "s"  domain  in  lieu  of  the  time  domain  is 

preferred  since  the  functional  form  of  typical  transfer  functions 

for  many  complex  processes  involving  heat  and  mass  transfer  have 

been  investigated  and  reported  in  the  literature.   See  for  instance 

references  [13],  [18]  and  [19].   An  alternative  route  for  getting 

an  approximate  transfer  function  is  to  plot  the  normalized  model 

response  on  semilog  paper.   The  equation  of  the  best  asymptote  to 

this  curve  yields  the  first  terra  in  the  approximate  equation  of  the 

response,  i.e.,  X(t)  =  a  e    *  +  a  e    3  + Next  the 

difference  between  this  asymptote  and  the  actual  response  is  plotted 

on  semilog  paper;  the  equation  of  the  best  asymptote  gives  the  second 

term  of  the  equation  for  X(t).   In  a  like  manner,  the  differences 

are  plotted  until  as  many  terms  as  needed  are  obtained  for  the  degree 


43 


of  accuracy  desired.   A  least  squares  fit  could  then  be  used  as  be 

fore 

X(s) 


fore  to  optimize  the  parameters  in  X(s)        with  respect  to 


observed' 

Transfer  Function  for  Top  Product  Response  to  a  Step  in  Feed 

The  model  contained  an  accumulator,  delays  in  condenser  and  re- 
boiler,  and  liquid  level  controls  in  reboiler  and  accumulator. 
The  transfer  function  G(s)  for  the  situation 


STEP  CHANGE     C, (s) 

IN  FEED    * ■ 


G(s) 


4V  J       ..  TOP  PRODUCT 


was  the  case  studied.   G(s)  was  assumed  to  be  of  the  form 


~rAS 
G(s)  = 


(T]Ls+1)(t2s+1) 


and  the  values  of  K,  t,»  t,  and  t9  were  estimated  graphically  as 
described  by  Moczek  et  al.  [13]  from  the  open  loop  time  responses. 
The  results  obtained  were  as  follows: 

K        Td  Tl  T2 

1       .3  tau      36.7  tau      4.1  tau 
or  units        units        units 

1       .08  min      9.8  min       1.09  min 

If  one  denotes  F(s)  as  the  perturbation  and  X(s)  as  the  response, 
the  G(s)  =  X(s)/F(s).   Thus,  for  a  step  in  C,  composition  of  the 
feed  of  0.05,  F(s)  =  .05/s  and  the  functional  form  of  X(s)  to  be 
used  with  equation  (24)  becomes 

x«>--    •05e"'3s 


i(36.7s+l)(4.1s+l) 


A  factor  of  100  was  used  to  expand  the  useful  X(t)  interval  so 
that  the  actual  form  of  X(s)  used  in  the  least  squares  fit  was 
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X(s/100) 

100       s  =  1,2, N 


The  least  squares  program  used  to  optimize  the  three  parameters 
(K,  t,  >  t9)  is  given  in  Appendix  V.   The  results  for  a  three  para- 
meter fit  were  K  =  .0536,  t.  =  36.2004  and  t  =  4.4585  or 


.      1.072  e",3s 

GQs;  "  (36.2004s+l)(4.4585s+l)  ^; 


Using  equation  (25),  X(s)  =  ,05/sG(s)  was  inverted  back  to  the  time 
domain  and  plotted  along  with  the  model  response  for  X(t).   The 
results  are  shown  in  Fig.  14.   The  fit  is  excellent  with  only  small 
deviations  at  large  values  of  elapsed  time--well  beyond  the  control 
region. 
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6.   Summary  and  Recommendations. 

A  dynamic  model  has  been  developed  which  performs  exceedingly 
well  under  the  assumptions  made  in  the  mathematical  development. 
Additional  effects  can  easily  be  incorporated  into  the  present  model 
once  the  mathematics  are  specified.   Preliminary  control  studies 
using  the  model  for  error  generation  indicate  single  loop  and  multi- 
loop studies  can  be  done  using  a  reasonable  amount  of  computer  time 
and  appropriate  transfer  functions  can  be  obtained  from  the  model 
responses.   All  computation  time  can  be  substantially  reduced  by 
converting  the  present  Fortran  60  programming  to  Fortran  63  or 
Fortran  4.   Thus,  the  model  can  be  a  useful  tool  for  design  and  con- 
trol work. 

Future  work  might  be  directed  toward  inclusion  of  pressure 
effects  in  the  column,  downcomer  effects,  reboiler  dynamics  and  a 
partial  condenser  to  extend  the  range  of  control  variables  and 
interactions  which  affect  the  system.   Data  on  an  actual  operating 
column  would  be  quite  helpful  in  this  respect,  so  that  the  model 
could  be  tailored  to  the  actual  plant  conditions.   In  this  way  control 
studies  and  model  responses  could  be  correlated  and  checked  against 
an  actual  operational  plant  in  hopes  of  improving  the  design  and 
control  loops. 
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APPENDIX   I 


DEVELOPMENT  OF    HYDRAULIC   EQUATIONS 


The  desired  functional  form  of  the  expression  relating  liquid 
holdup  to  column  flows  is  W   =  f(Lj).   The  following  development 
is  adapted  from  an  analysis  for  perforated  plates  done  by  Peiser  and 
Grover  [1], 


"»-; 


% 


LJ 


AT; 


Yi'-, 


M*  I 


\ 


AT3-= 
h  = 

1.  H 

J 

g  * 

V  = 
L  = 


weir  height (ft) 
=  liquid  height  on  the  tray(ft) 

=   liquid   height   in  downcomer(f t) 

2 
=  cross  sectional  area  of  downcomer(ft  ) 

2 
=  cross  sectional  area  of  tray(ft  ) 

3 
=  liquid  molar  density(raoles/f t  ) 

=  height  of  the  weir  crest(ft) 

width  of  the  weir(ft) 

2 
gravitational  constant(ft/sec  ) 

molar  vapor  f low(moles/time) 

molar  liquid  f low(moles/time) 
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c  =   correlation  constant  in  the  Francis  weir  equation 
cf)   =  liquid  flow  resistance  in  downcomer 
en  .  =  liquid  flow  reversal  resistance  in  downcomer 
cp   s  pressure  drop  for  vapor  flow  and  liquid  head  on  tray  abobe 
U   =  liquid  velocity  in  downcomer 
C,  .  =   locally  defined  constant 
A^  s  area  under  downcomer 

v.  =  vapor  velocity  at  tray  (j  +  1) 

A  =  perforation  area 

N  =  total  number  of  perforations  at  tray  (j  +  1) 

C  =  orifice  discharging  coefficient  (approx.  0.6) 

M  =  molecular  weight  of  liquid 

M  =  molecular  weight  of  vapor 

Using  the  Francis  weir  formula  for  large  straight  weirs 

Lj  =  kh  (1) 

where 


k  =  *m  ijPLj 


From  the  diagram 

h  =  cpTj  -  %i  (2) 

Using  a  pressure  balance 

%r  \r  14Kt-!{- cuLi+°  <4) 
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where 


1-5    /  a         ^ 


'Lj        .36(2g)^Lj     VDj' 


In  a  similar  manner 


Li! pvi    "vi 

fPJ   "  CDj+1(2g)        PLj     ^.  «!j+l 


(5) 


vJ=4j1VpJ 


j+1 


(6) 


Therefore, 


*Pj  =  CVj+lVj  +  ^Tj+l 


(7) 


where 


CVj+l 


**V    I       I   (C  AN  )2  , 

L28MLPLPvJj/ 


Finally 


*D2   '  *Ej  +  CLjLj+l  +  Vj     +  «rjrt 


n 


J    "    %j     !"    k   Lj 


* 


(3) 
(9) 


where  k1   =  (I-y$ 


Liquid  holdup3 


WLi   =   <PLATVl  +   (plYV 


(10) 


3  3 


substituting  equations  (8)  and  (9)  into  equation  (10)  yields 

(ID 


j  =  pLjATjkj  +  A|j  +  PLAj  [«Ej  +  *Tj+l  +  CLjLj+l  +  CVj+lVj2] 


This  expression  is  for  a  perforated  plate.   To  simplify  the 
analysis  for  a  bubble  cap  tray,  it  will  be  assumed  that 

(a)  volumetric  content  of  hydraulic  gradient  and  downcomer 
equa{b  volumetric  content  of  risers  and  caps. 

(b)  unidirectional  flow  predominates. 

(c)  Francis  weir  equation  holds. 

Let  cp_j  =   effective  liquid  depth  on  plate  j  and  A  =  effective 
surface  area  over  which  cp   applies.   Then 

WLj  =  ^j^Tj  (12) 

Substituting  equations  (1)  and  (2)  into  (12)  yields 

wlj  -  ^lja[%  +  <yk>  'J  (") 

Converting  all  flows  to  flow/unit  feed  flow  and  using  a  holdup 
ratio,  R^  =  wLj/(wLr)ss»  equation  (13)  becomes 


where 

A      -ft,  „  _    A 


3a  <Mss  yp^H  J 


Lr^SS  ^WLr' 


Differentiating  equation  (14)  yields 
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APPENDIX  II 

FORTRAN  60  PROGRAMS  USED  TO  SIMULATE  COLUMN  RESPONSES  ON  A 
CDC- 1604  DIGITAL  COMPUTER 


Nomenclature  pertinent  to  all  programs: 


X 
XIF 


YIF 


XISSF 


YISSF 


single  component  mole  fraction  in  stage  liquid 

single  component  mole  fraction  in  entering 
liquid  feed  (FL> 

single  component  mole  fraction  in  entering 
vapor  feed  (F  ) 

single  component  mole  fraction  in  liquid  side 
stream  feed  (Lgs  ) 

single  component  mole  fraction  in  vapor  side 
stream  feed  (Vssp) 


TEMP  temperature  of  stage  liquid 

VAP  vapor  flow 

QUID  liquid  flow 

FL  liquid  feed  entering  a  stage 

FV  vapor  feed  entering  a  stage 

SSFL  liquid  side  stream  product  leaving  a  stage 

SSPV  vapor  side  stream  product  leaving  a  stage 

VAPMOH  vapor  molar  enthalpy 

QUIMOH  liquid  molar  enthalpy 

HFL,  HFV,  HSSFL,  HSSFV  enthalpies  of  streams  indicated 

RW  liquid  holdup  on  a  stage 

RWDOT  liquid  holdup  derivative 


ENTHK,  ENTHL 


ENTHV,  ENTHW 


constants  in  the  liquid  enthalpy  equation 

(h.  =  a.T.  +  b.) 
i    i  J    i 

constants  in  the  vapor  enthalpy  equation 

(H.  -  aJT.  +  b!) 
i    i  J    i 
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A,  B,  C 

ALPHA,  BETA,  GAMMA 

DELTAU 

BPERR 
DELERR 

RHOQD 

QC,  QR 

NC,  NS,  JFEED,  NCNTROL 

WEIR 

HW 

AREA 

FD 

CT 

GR 

RHOL 

HT 

WLRSS 

ALPH,  BET 

DERIVM 

DERIVH 

QUNBAL 

NTALLY 


constants  in  the  equilibrium  equation 
(K.  =  exp(Ai/Tj  +  B^^  +  C^) 

constants  in  the  Runge-Kutta-Gill  integration 

computation  interval  for  integration 

error  limit  allowed  in  bubble  point  calculation 

error  limit  allowed  in  convergence  test 
h\eqn(8)|  <e) 


(|h.eqn(20) 


individual  component  densities,  p 
condenser  and  reboiler  duties 


Li 


number  of  components,  number  of  stages,  feed 
stage,  number  of  increments  over  which  inte- 
gration is  to  proceed 

weir  length 

weir  height 

effective  tray  area 

feed  rate(F) 

constant,  c,  in  the  Francis  weir  equation 

acceleration  of  gravity,  g. 

total  liquid  density 

effective  height  of  liquid  on  a  tray 

<WLr>SS 

constants  a  and  8  in  equation  (9) 

Ui 

h.  from  equation  (20) 

•  J 

h .  from  equation  (8) 

number  of  At  intervals  elapsed 


The  following  programs  are  for  the  model  with  varying  boundary 
conditions.   All  programs  utilize  a  total  condenser  and  the  hydraulic 
interactions,  equations  (9). 
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C  THIS  PROGRAM  HAS  AN  ACCUMULATOR*  BUT  NO  CON.  OR  REBOIL.  DELAYS 


PROGRAM  DYNAMIC 

DIMENSION  X(10.50).TX(10.50)»XIF(10.50).YIF< 10. 50 ) .X  I SSF ( 10.^0) ,Y I 
1SSF( 10.50).ENTHK( 10) .EKC( 10) . 

2    ENTHL( lO).ENTHV(lO)  .ENTHW(IO) .EM  10) .ALPHA ( 3 ) .BE T A ( 3 ) .GAMMA ( 3) . 
3DERIVM( 10) »TEMP(50).VAP<50),QUID(50)»FL(50).FV(50).SSFL(50) .SSFV(5 
40).SSPL(50) »SSPV(50) .SUM( 50 ) »VAPMOH( 50 ) .QUI MOH ( 50 ) .DER I VH( 5  0 ) .QUNB 
5AL(50) »HFL( 50).HFV(50) .HSSFL(50) .HSSFV(50) »RW( 50) »RWDOT(50) , 
6RHOOD( 10) .CON( 10) .FEED (10) .QUI DX ( 1 0 ) ♦ A ( 10 ) ♦ B ( 10 ) . C < 10 ) .DLYX < 10 ) 

COMMON  QUIDX.A.B.C 

EQUILKF(A.B.C»T)=EXPF(A/(T+460.0)+B+C«(T+460.0>) 

THALPF(Y.Z.T)=Y*T+Z 

READ  2.NC.NS.JFEED.NCNTROL 
2  FORMAT (41 5) 

READ  333. ( ALPHA ( I ).BETA(I) .GAMMA (I ).I=1.3) 
333  FORMAT(6F12.6) 
1  FORMAT(8F10.4) 

JLIM=NS+2 


READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ  1 

READ  1 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ1 


READ 
READ 
READ 
READ 
READ 
READ 
READ 
READ 
READ 


UX(I.J).I=1.NC).J=1.JLIM) 

<(XIF( I.J) »I=1.NC).J=2.NS) 

( (YIF( I,J).I=1»NC) »J=2.NS) 

( (XI SSF ( I .J). 1  =  1. NO  »J=2.NS) 

( (YISSFf  I,J).I«1.NC) ,J=2.NS) 

(TEMP( J) .J=1»JLIM) 

(VAP(J),J=1.JLIM) 

(QUID(J) .J-l.JLIM) 

(FL( J) .J=2.NS) 

(FV( J) »J=2.NS) 

(SSFL(J) .J=2.NS) 

(SSFV(J) .J=2.NS> 

(SSPL(J) .J=2.NS) 

(SSPV( J) .J=2.NS) 

(VAPMOH( J).J=1.JLIM) 

(QUIMOH( J) # J>ltJLIM) 

(HFL( J).J=2.NS) 

(HFV( J).J=2.NS) 

(HSSFL( J).J=2.NS) 

(HSSFV( J).J=2»NS) 

(RW( J) »J=2»JLIM) 

RWDOT( J) .J=2.NS) 

(ENTHM I ) »I=1»NC) 

(ENTHL( I).I=1.NC) 

(ENTHV( I ).I=1.NC) 

(ENTHW( I >.I=1.NC) 

(A( I ).I=1,NC) 

<B(  I  )  .1  =  1. NO 

(C(  I  )»I  =  l.NO 

DELTAU.BPERR.DELERR 

(RHOQD( I ) .1=1. NO 


OR  INITIAL  =  13687.026 

QR  UNDER  PERTURBED  CONDITIONS  IS  NOW  SET  AT  13247.239 


OR  =  13247.23948193 
WEIR=4.62 
HW=0.25 
AREA=21.49 
FD=0. 27626 
CT=0.42 
GR«32.2 
RHOL=0.0 
DO  100  1=1 ,NC 
100  RHOL=RHOL  +  X ( I  »  JFEED) «RHOOD( I ) 
C0NST=CT*SQRTF(2.*GR)*WEIR 
DENOM=CONST*RHOL 

HT=HW  +  (FD/DEN0M)**.6667*(QUID( JFEED) )**. 6667 
WLRSS=RHOL*AREA*HT 
ALPH=AREA*HW/WLRSS 

BET=AREA/WLRSS*(FD/C0NST)**.666  7 
LABEL  =8HALPH  BET 
PRINT  6.LABEL.ALPH.BET 
6  F0RMAT(A8*9F12.6) 
DO  9  J=3»NS 
RHOL  »0.0 
DO  13  1=1. NC 
13  RHOL=RHOL  +  X (  I  ,  J ) *RHOOD(  I ) 
9  RW(J)=ALPH*RHOL  +  BET* ( RHOL**. 3333 ) * (QU I D( J ) **.6667 ) 
PRINT  1000»<RW( J) »J»1»JLIM) 
1000  F0RMATU5H  HOLDUP  RATIOS  ./»12F10.6) 
PHI  =  0.1*DELTAU 
NTALLY=0 
NPRINT=0 
NPRINT  =  0 

C  FSN  445  SEQUENCE  IS  FOR  STAGES  BELOW  FEEDSTAGE*  444SEQUENCE  IS  ABOVE 


4  JPATH=1 

JVARY=JFEED 

JLOW=JFEED 

JHIGH=NS 

NUM=+1 

GO  TO  444 
445  JLOW=3 

JHIGH=JFEED-1 

JVARY=JHIGH 

NUM=-1 
444  DO  26  M=JLOW.JHIGH 

J=JVARY 

JVARY=JVARY+NUM 

LTIMES=1 


C  THIS  SECTION  THRU  FSN  12  IS  R-K-G  INTEGRATION  FOR  A  GEN.  STAGE 
44  DO  11  I=1»NC 


TX(  I.J-1)=X(  I  .J-l) 
TX( I»J)=X( I »J) 

11  TX( I .J+l) =X( I t J+l ) 
DO  12  1=1. NC 
QUE=0.0 

EQK  =  EOUILKF(A( I  ) .S(  I  ) .C(  I  )  ♦  TEMPI  J)  ) 

TERM=VAP( J)*EGK+QUID( J ) 

CON( I )=(VAP( J-l )*EQUILKF( A( I ) , 3< I ) » C ( I ) • TEMP ( J-l ) ) *TX ( I • J- 1 ) )+QUID 
1  (J+l )«TX(  I , J+l ) 

FEED(I)=FL(J)*XIF(  I ♦ J ) +FV ( J ) *Y I F (  I ♦ J )  +  SSFL  (  J ) «X  •  SSF  (  I • J ) +SSF V( J ) *Y 
1ISSF( I.J) 

OUTPUT =TX(  I  ,J)*(SSPL( J  )+SSPV ( J ) *EQK ) 

CAY  =  (  CON(  I  )  -TX<  I  .  J)*f  TERM+RWDOK  J)  )  +FEEDI  I  )  -OUTPUT  )*DELTAU/RW(J) 

DO  5  K=1.3 

TX( I»J)=TX(  I  ♦J)+GAMMA(K)*(CAY-QUE) 

QUE=ALPHA(K  )*CAY+BETA(K)*QUF 
5  CAY=(CON(  I  >-TX(  I . J )* { TERM+RWDOT ( J  )  )  +FEED(I) -OUTPUT ) *DELTAU/RW( J) 

12  TX( I»J)=TX( I ♦J)+CAY/o.0~QUE/3.0 
31  5UM(J)=0.0 

DO  7  !=1.NC 

7  5UM( J)=5UM{ J)+TX( I ,J) 
DO  8  I=1»NC 

TX( I »J)=TX( I »J)/SUM( J) 

8  OUIDX( I )=TX( I.J) 

CALL  BUBPT(TEMP(J) .NC.BPERR) 

VAPMOHM  J)=0.0 

0UIM0H( J)=0,0 

DO  15  1=1. NC 

EK( I )=EQUILKF(A( I ) »B( I ) .C< I ) »TFMP( J) ) 

DERIVM( I )  =  (CON(  I )-TX(  I  .J)*(EK( I ) *VAP ( J ) +QU ID ( J > +RWDOT ( J ) )+FEED(  I)- 
1  TX(I.J)*(SSPL( J)+SSPV(J)*EK( I ) ) )/RW( J) 

QUIMOhM  J)=QUIMOH( J ) +  THALPF ( ENTHK ( I  ) .ENTHL( I ) .TEMPI  J)  )*TX(  I.J) 
15  VAPMOH<  J)=VAPMOH(  J  ) +THALPF  (  ENTHV  (  I  ),ENTHW(  I  ).TEMP(J)  )*TX(  I  »J)«EIC(  I 
1) 

TERM1=0.0 

TERM2=0.0 

TERM3=0.0 

TERM4=0.0 

DO  21  1=1. NC 

TERMl=TERMl+( ENTHK ( I ) *TEMP ( J ) +FNTHL ( I ) )*TERIVM( I ) 

TERM2=TERM2+EK( I )*DERIVM( I ) 

TERM3  =  TERM3+EK(  I)*TX(  I.J)*(A<I)/(  ( TEMP ( J ) +460 .0 ) **2 ) -C (  I )  ) 
21  TERM4=TERM4+ENTHK(  I )*TX(  I.J) 

DERIVH{ J)=(TERM1+(TERM2*TERM4)/TERM3) 

FDTHALP=VAP( J-l ) *VAPMOH ( J-l ) +OUID < J+l ) *QUI MOH ( J+l ) +FL ( J ) *HFL ( J ) +F V 
1 ( J)*HFV( JJ+SSFL ( J)*HSSFL( J)+5SFV( J)*HSSFV( J) 

GO  TO  ( 33,33.221) .JPATH 
3  3  OUNBAL( J)=( FDTHALP-QU I MOH ( J ) * ( SSPL ( J ) +RWDOT ( J ) +QU ID ( J ) )-VAPMOH( J)* 
1 (SSPVt J)+VAP( J) ) )/RW( J) 


C   THE  NEXT  STATEMENT   COMPUTES  THE  DERIVATIVE  DIFF  FOR  CONVERG.  TEST 


DIFF=DERIVH( J)-QUNBAL( J) 
IF(ABSF(DIFF)-DELERR)2  5.25.2  2 


C  FSN  22    THRU  26   IS  ADJUSTMENT  OF  COEFF.  BY  REGULI  FALSI 


22  LTIMES=LTIMES-1 
RHOL=0.0 

DO  101  I=1.NC 
101  RHOL=RHOL+TX( I»J)*RHOQD( I ) 
EXP=RHOL**.3333 
IF(LTIMES)  24,23.23 

23  RWO=RW(J) 
RWREF=RW( J) 

RW<J)  =  RW(J)  +  PHI 
RWDOT( J)=.l 

OUID(J)=( (RW(J)-ALPH*RHOL)/(BET*EXP) )**1.5 
SUMFDS=VAP( J-1)+QUID{ J+1)+FL( J)+FV( J)+SSFL( J)+SSFV( J) 
VAP( J)=SUMFDS-QUID(J)-RWDOT( J) 
35  DIFFO=DIFF 
GO  TO  44 

24  SLOPE=(DIFF-DIFFO)/(RW( J)-RWO) 
RWNEW=RW( J) -D IFF/SLOPE 
RWDOT( J)=(RWNEW-RWREF)/DELTAU 

QUID(J)=< (RWNEW-ALPH*RHOL)/(BET*EXP) ) ** 1 • 5 
VAP ( J ) =SUMFDS-QU ID ( J ) -RWDOT ( J ) 
DIFFO=DIFF 
RWO=RW< J) 
RW( J)=RWNEW 

IF(LTIMES+30)  37,44.44 
37  LABEL=8HEND  TEST 

PRINT  6,LABEL,J,QUNBAL( J) ,DERIVH( J) ,OU I D ( J ) , VAP ( J ) ,RW( J ) , RWDOT ( J ) 
GO  TO  30 

25  DO  26  1=1, NC 

26  X(  I,J)=TX( I ,J) 
JPATH=JPATH+1 

GO  T0<199, 199,201) ,JPATH 


C  FSN  199  THRU  219  INTEGRATES  THE  CONDENSER/ACCUM  SECTION 


199  J=NS 

DLYV=VAP( J) 

DLYHV=VAPMOH( J) 

DLYHL=0.0 

DO  502  1=1. NC 

EK( I )=EQUILKF(A( I ),B( I ),C( I ),TEMP( J) ) 

DLYXd  )=EK(  I  )*X(  I, J) 
502  QUIDX( I )=DLYX(  I  ) 

TEMPCON=TEMP( J) 

CALL  BUBPT(TEMPCON,NC,BPERR) 

DO  504  1=1, NC 
50  4  DLYHL=DLYHL+THALPF(ENTHK( I ),ENTHL( I ) ,TEMPCON ) *DL YX ( I ) 

CONDMOH=DLYHL 

QC=DLYV    *(DLYHV    -DLYHL    ) 

C0NLI0=0UID( J+l )+QUID( J+2) 


LTIMES»1 

RWREF=RW< J+l) 

RWDOT< J+l )=DLYV    -CONLIQ 

RW( J+l)=RWREF+RWDOT( J+l )*DELTAU 

205  DO  206  1=1. NC 

206  TX( I.J+1)=X< I. J+l) 
DO  207  1=1. NC 
QUE=0.0 

CAY=(DLYV*DLYX( I  )-( CONL IQ+RWDOT < J+l ) )*TX( I . J+l ) )*DELTAU/RW ( J+l ) 

DO  208  K=1.3 

TX< I»J+1)=TX( I.J+1 )+GAMMA(K)*(CAY-QUE) 

QUE=ALPHA(K)*CAY+BETA(K)*QUE 

208  CAY=<DLYV*DLYX( I ) - ( CONL IQ+RWDOT ( J+l ) )*TX< I. J+l) )  *DELTAU/RW ( J  +  l ) 

207  TX(  I.J+1)=TX(  I.J+D+CAY/6.0  -  QUE/3.0 
SUM( J+1)=0.0 

DO  209  I=1.NC 

209  SUM(J+1)=SUM(J+1)+TX< I.J+1) 
DO  210  1=1, NC 

TX( I.J+1)=TX( I.J+1) /SUM (J+l) 

210  OUIDXI I)=TX( I.J+1) 

CALL  BUBPT(TEMP(J+1) .NC.BPERR) 

QUIM0H( J+l) =0.0 

DO  211  I-ltNC 

EKC(I)=EQUILKF(A(  I  ).B(  I).C(I  ).  TEMP  (J  +  l)  ) 

DERIVM( I)=(DLYV*DLYX( I )-( CONL IQ+RWDOT ( J+l ) )*TX( I » J+ 1  )  ) /RW( J+l ) 

211  QUIMOH( J+l)=QUIMOH( J+l ) +THALPF ( ENTHK. ( I ) »ENTHL( I ) » TEMPt J+l ) ) *TX< I . J 
1  +  1) 

TERM1=0.0 

TERM2=0.0 

TERM3=0.0 

TERM4=0.0 

DO  212  1=1, NC 

TERM 1= TERM 1  +( ENTHK ( I ) *TEMP ( J+l )+ENTHL ( I ))*DERIVM( I ) 

TERM2=TtRM2  +EKC ( I ) *DER I VM ( I ) 

TERM3=TERM3+EKC( I ) *TX( I . J+l ) *( A( I ) / ( ( TEMP( J+l ) +460. ) **2 ) -C ( I)) 

212  TERM4=TERM4+ENTHK( I )*TX( I.J+1) 

DERI VH( J+l) =(TERM1+(TERM2*TERM4)/TERM3) 

QUNBALI J+1)=(DLYV    *C0NDM0H- ( CONL IQ+RWDOT ( J+l ) )*QUIM0H( J+l ) )/RW( J 
1  +  1) 
216  DO  219  1=1. NC 
219  X( I.J+1)*TX(I.J+1) 

GO  TO  445 


C  FSN  201  THRU  228  INTEGRATES  THE  RFBOILER  SECTION 


201  J=2 

LTIMES=1 
GO  TO  44 

221  QUNBAL ( J ) = ( FDTHALP-V AP ( J ) *VAPMOH ( J ) -QU I MOH( J ) * ( QU I D( J ) +RWDOT ( J ) ) +Q 
1R)/RW( J) 

DIFF=DERIVH( J) -QUNBALI J) 

IF ( ABSF ( D IFF )-DELERR) 22 5, 22 5.2 22 

222  LTIMES=LTIMES-1 
IF(LTIMES)224, 223,223 


223  RWO=RW(J) 

RWREF=RW< J) 

RW( J)=RW< J)+0.01*DELTAU 

RWDOT( J)=0.01 

VAP( J)=QUID( J+1)-QUID< J)-RWDOT( J) 

DIFFO=DIFF 

GO  TO  44 
2  24  SLOPE=(DIFF-DIFFO)/(RW( J)-RWO) 

RWNEW=RW( J) -D IFF/SLOPE 

RWDOT( J)=(RWNEW-RWREF)/DELTAU 

VAP( J)=QUID( J+1)-QUID( J)-RWDOT( J) 

DIFFO=DIFF 

RWO=RW( J) 

RW( J)=RWNEW 

IF (LTIMES+30) 2  26*44,44 

226  PRINT  227,QUNBAL(J),DERIVH(J) 

227  FORMAT( 13HREBOILER  LOOP . 2F20 .8 ) 
GO  TO  30 

225  DO  228  1=1. NC 

X( I.J-1)=TX( I, J) 

228  X( I ,J)=TX< I ,J) 


C  FSN253  THRU  STOP  IS  THE  DATA  PRINTOUT  SECTION 


253  NTALLY=NTALLY+1 

IF(NTALLY-NPRINT)4,17,17 

17  PRINT  18,NTALLY.DELTAU 

18  FORMAT( 17H1INCREMENT  NO.  =  I 5/34H0 I NDEPENDENT  VARIABLE  INCREMENT  = 
1  F10.4) 

PRINT  19 

19  FORMAT( 116H0STG    MOL  FRA  SUM    H  UNBALANCE   H  DERIVATIVE    TEMPER 
1ATURE    LIQUID  FLOW   LIQ  ENTHALPY    VAPOR  FLOW   VAP  ENTHALPY      ) 

JSTG=0 
JLIM=NS+1 
DO  27  J=2»JLIM 

PRINT  20,JSTG.SUM( J) ,QUNBAL( J) ,DERIVH( J) . TEMP ( J ) , OU I D( J ) .QUI MOH ( J ) 
1 »VAP< J) .VAPMOHt J) 

20  FORMAT( I4.8E14.8) 

27  JSTG=JSTG+1 
PRINT  271.QC.QR 

271  FORMAT ( 18H0CONDENSER  DUTY  =  El5. 8/ 17H0REBOI LER  DUTY  =  E15.8) 
JSTG=0 

DO  29  J=2.JLIM 
PRINT  28.JSTG.(X( I ♦ J ) . I = 1 .NC ) . RW( J ) »RWDOT( J) 

28  FORMATC9H0       STAGE  NO.  =  I  3/  (  1 0F13.  10  )  ) 

29  JSTG=JSTG+1 
NPRINT  =  NPRINT  +  1 
IF(NTALLY-NCNTROL)  4,30.30 

30  STOP  77777 
END 

C        SUBROUTINE  FOR  BUBBLE  POINTS 
C 

SUBROUTINE  BUBPT ( T ,L , BPERR  ) 

DIMENSION  GUIDX(10),A(10) ,B(10) ,C(10) 


COMMON  OUIDX.AtB.C 

EQUILKF<A.B.C.T)=EXPF(A/(T+460.0>+B+C*<T+460.0)  ) 
KTIMES=1 

3  SUMY=0.0 
DO  4  I=1,L 

Y=EQUILKF(A(I).B< I  >»C<  I )»T)*Q(JIDX< 1 ) 

4  SUMY=SUMY+Y 
IF(ABSF(SUMY-1.0)-BPERR)8.8t3 

5  KT1MES=KTIMES-1 
IF(KTIMES)  7,6.6 

6  SUMYO=SUMY 
TO=T 

T=T+10.0 
60  TO  3 

7  SLOPE=(SUMY-SUMY0)/(T-T0) 
TN=< (1.0-SUMY) /SLOPE )+T 
SUMYO=SUMY 

TO=T 
T=TN 
GO  TO  3 

8  RETURN 
END 
END 


1.3 


PROGRAM  DYNAMIC 
THIS  PROGRAM  HAS  NO  DECAYS  OR  ACCUMULATOR 

DIMENSION  X<10,50)  ,TX (10,50 ),XIF( 10,50 ),YIF( 10,50 ),XiSSF(lQ,50)#YJ 
lSSF<10,5o ) ,EN7HK<10 ) - 

2    ENTHLC10) .EN7HV{10),EN7HW(10) .EK(10)# ALPHA(3).BETA(3)#GAMMA<3). 
3DERIVM(io)#TEMP(5o ) ,VAP(5o)»QUID(5o),FLC50)»FV(5o),SSFL(5o),SSFV(5 
40)*SSPL<50),SSPV(5U),SJM(50),VAPMOH<50>.QUIMOH(50) . DER I VH ( 5 0 ) , QUNB 
5AL(5o>,HFL(5o),HFV(5o>,HSSFL(5o),HSSFV(5o),Rrf(5o)»RWDOT(5o). 
6RHOOD(10),CON(10),FEED(10),GUIDX(10),A(10)#B(10),C(10) 

COMMON  Q  U  I  D X . A , B , C 

EOUILKF(AJB,C,T)=EXPF(A/(T+460.0)+B+C-(T+460.0)) 

THALPF(Y,Z, T)=Y*T+Z 

READ  2,  NO  NS,  JFEED,NCMROL 
2  FORMAT (4 15) 

READ  333» (ALPHA ( I), 8E7AC ), GAMMA ( I), 1=1,3) 
333  FORMAT (6F12. 6) 
1  FORMAT (3F10 .4) 

JLIM  =  NS-r2 

READ  l,((X(I,J),I=l,NO,J  =  l,JLIM) 

READ  1,  (  (XIF(  I,  J),  I  =  l/-NO,J  =  2,NS) 

READ  1,  (  (YIF(  I,  J),  I=l,NO,J  =  2*NS) 

READ  i,  (  (XISSFU,  J),  I=l,NO#  J  =  2/NS) 

READ  1,  (  (YISSF(  I,J>,  I=l,NO/ J  =  2»NS) 

READ  1, ( TEMP( J) , J  =  l, JLIM) 

READ  1, (VAP( J) , J=l, JLIM) 

READ  1, (QUID( J) , J=l, JLIM) 

READ  1,  (FL( J) ,J=2,NS) 

READ  1,  <FV( J) ,J=2,NS) 

READ  1, (SSFL( J) , J=2,NS) 

READ  1, (SSFV( J) , J=2,NS> 

READ  1, (SSPlC J) , J=2,NS) 

READ  1,  (SSPV( J) , J=2,NS) 

READ  1, ( VAPMOHC J) , J=l, JLIM) 

READ  1, (QUIMOHC J) , J=l, JLIM) 

READ  1, (HFL( J) , J=2,NS) 

^EAD  1, (HFV< J), Js2# NS> 

^EAD  1, (HSSFL( J) , J=2,NS) 

READ  1,  (HSSFV( J)  ,  J=2,NS) 

READ  1, (RW( J) , J=2» JLIM) 

READl. (RWOOT( J) , J=2,NS) 

READ  1,  (ENTHK(  I  ),  1=1,  NO 

READ  1/ (ENTHLCI). 1=1, NO 

READ  1,(EN1HV(I),I=1,NC) 

READ  1,  (ENTHW<  I).  1=1, NO 

READ  l,(A(I),I=i,NO 

READ  1, (B(I ) , 1=1, NO 

3EAD  1,(C(I), 1=1,  NO 

READ  1,DEL.TAU,BPERR,DELERR 

3EAD  1,(RHOQD(I),I=1.NO 


QC  ^    /3Z47.  239*0/93 


,;e:R  =  ^. 62 

r;  'fi  =  C  .  2  5 

AREA=2i.49 
r~D  =  0.  27626 


CT=0.s2 

GR=32-2 

RHOl=0.0 

DO  10  0  1=1, NC 
100  RHOL=RHOL  ♦  X ( I , JFEED ) *RHOQD C I  ) 

CONST =CT«SGRTF( 2. *GR)« WEIR 

DENOM=CONST*RHOL 

HT  =  HW  *  iFD/DEN0M)**.6667-(QUlD< JFEED)  )**.6667 

WI_RSS  =  RHOL*AREA*HT 

ALPH=AREA*HW/wLRSS 

BET  =  AREA/WLRSS*(FD/C0N,ST)**,6667 

LABEL  =8hALPH  BET 

PRINT  6, LABEL, ALPHABET 
6  F0RMAT(A3,9Fi2.6) 

DO  9  J=3,NS 

RHOL  =0-0 

DO  13  1=1, NC 
13  RHOL=RHOL  +  X ( I , J ) *RHCQD (  I  ) 
9  RW( J)=ALPH*RHOL  ♦  BET* ( RHOL**  .  3333  )*( QUI D < J >**. 6667 ) 

PRINT  1QC0, (RW( J), J=l, JLIM) 
1000  FORMATdSH  HOLDUP  RATIOS  ,/,12F10.6) 

PHI  =  0.1*DELTAU 

NTALLY=0 

NPRINT=0 

4  JPATH=1 
JVARY=JFEED 
JLOW=JFEED 
JHIGH=NS 
NUM=+1 

GO  TO  444 
445  JL0kl  =  3 

JHIGH=JFEED-1 

JVARY=JHIGH 

NUM=-1 
444  DO  26  MsJLOW, JHIGH, 

J=JVARY 

JVARY=JVARY+NUM 

LTIMES=1 
44  DO  11  1=1, NC 

TX(  I, J-1)=X(  I , J-l) 

TX( I. J)=X( I, J) 

11  TX( I, J+l)sX( I,J+1) 
DO  12  1=1, NC 
QUE=0-0 

EQK=EQUILKF(A(I),B(I),C(I),TEMP(J>) 
TERM=VAP( J)*EQK*QUID( J) 

CON( I )=(VAP( J-1)*E0UILKF(A( I ),8( I ),C( I ), TEMP ( J-l )  )*TX ( I, J-l) )*QUID 
1<J-1)»TX< I, J  +  l) 

FEED(I)=FL(J)*XIF(I,J)+FV(J)*YIF(I,J)*SSFL(J)*XISSF(I,J)*SSFV(J)*Y 
1ISSFCJ) 

OUTPUT  =  TX(I,J)*(SSPL<J)«'SSPV<J)*EQK) 

CAY=(CON(I)-TX(I,j)*(TERM  +  RWDOT(J))  +  FEEDM  > -OUTPUT )* DEL TAU/RW ( J ) 

CO  5  K=i,3 

TX(  2.  J)sTX(  I  , J)*GAKMA(K>  -(CAY-QUE) 

CUE  =  ALPHA(K)*CAY  +  3t:~A(K)*QUE 

5  CAY=(CON(I)-TX( I, J)* i TEHM+RdDOT ( J > )  -FEED ( I  ) -OUTPUT ) *DELT AU/RW ( J ) 

12  VX(  I,  J)  =  TX<  I.  J)*CAY/6.'0-CUE/3.0 
31  SJM(J)=Q.O 

uO  7  i=l.NC 
7  SUM( J)=SJM( J)*TX( I , J> 
DO  8  I=1,NC 
TX(  I  ,  J)sTX(  I,  J)/SUM(  J)        f,5 


8  QU  I  DX  <  I  )=TX( I, J) 

CALL  3J8PT(TEMP<J),NC.BPERR> 

VAPMOHC J)=0.0 

QUIMOHC J)=0  .0 

DO  15  1=1, NC 

EKU)=EQUILKF(A<I  ),9(  I  ),C(  I  ),T£MP(J>) 

DERIVMd )s(CON(I >-TX( I, J)*(EK(i )*VAP( J)*QUID(J)*RWDOT(J) )*FEED<I)- 
1    TX(i.J)*(SSPL(J)+SSPV(J)vEK(;)))/Rw(J) 

QUIMOH<J)sQUIMOH(J>+THALPF(ENTriK{ I ),ENTHL(I >»TEMP(J))*TX<I#J) 
15  VAPMOH(J)=VAPMOH<J)+ThALPr"<ENTHV(I),ENTHW<I),TEMP(J))*TX(I#J)*EK(I 
1) 

TERM1=C.0 

TERM2=0  .0 

TER.vi3  =  0.0 

T£RM4=0.0 

DO  21  1=1, NC 

TERMi=TERMi+(£NTHK(I ) *TEMP ( J )+ENTHL ( I ) ) *DERI VM ( I ) 

TERM2  =  TERM2-'-EK(I)*DERIVM(l) 

T=RM3=TERM3+EK(I)*TX(I,J)*(A(I)/((TEMP(J)*460.0)**2)-C(I)) 

21  T£R,M4  =  TERM4  +  ENTHK(  I  )*TX(  I,  J) 
DERIVH(J)s(TERMl+(TERM2*TERM4)/TERM3) 

FDTHALP=VAP( J-l)*VAPMOH( J-1)*QUID( J*1)*QU1M0H( J*1)*FL< J)*HFL( J)*FV 
1(J)*HFV{J)+SSFL(J)*HSSFL(J)+SSFV(J)*HSSFV(J) 

GO  TO  (33,33,221), JPATH 
33  QUN3AL(J)=(FDTHALP-QUIM0H(J)*(SSPL(J)*RWD0T(J)*QUID(J) )-VAPMOH(J>* 
1(SS?V( J)+VAP( J) ) )/RW( J) 

DIFF=DERIVH( J)-QUN8Al< J) 

IF(A3SF(DIFF)-DELERR)25,25,22 

22  LTIMES=LTIMES-1 
RHOL=0 . 0 

DO  101  1=1, NC 
101  RHOL=RHOL+TX( I, j)*RHOGD( I ) 
EXP=RH0L**.3333 
IF(LTIMES)  24,23,23 

23  RwO=RW(J) 
RWREF=RW( J) 

RWCJ)  =  RW(J)  +  PHI 
RWDOT( J)  =  .l 

QUID(J)=((RW( J)-ALPH*RH0L)/<BET*EXP))**1.5 
SUMFDS=VAP(J-1)+QUID(J+1)+FL(J)+FV(J)+SSFL(J)+SSFV(J) 
VAP(J)=SUMFDS-QUID(J)-RWDOT( J) 
35  DIFF0=DIFF 
30  TO  44 

24  SLOPE=(DiFF-DIFFO)/(RW(J)-RWO) 
RWNEw=RW( J)-DiFF/SLOPE 
RrtDOTC J)=(RWNEW-RWREF)/DELTAU 

3UIDC J)=( (RWNEW-ALPH*RH0L)/(BET*EXP))**1.5 
VAP(J)=SuMFDS-QUlD(J)-RWDOT(J) 
D  i  F  F  O  =  D  I  F  F 
Ra'0  =  RW(  J) 
R  M  <  J  )  =  R  a  iM  E  W 

IFCLTIMES+30)  37,44,44 
37  lABEL=8HEND  TEST 

PRINT  6,  LABEL, J,QUNBAL( J) , DERIVH(J),QUID(J),VAP(j;,RW(J),RWDOT(J) 
GO  TO  30 

25  DO  26  1=1, NC 

26  X( I , J)=TX(  I  , J) 
J?ATrt  =  JPATti*l 

GO  T0(199, 199,201), JPATH 
199  J=NS 

QUlD(J+2)=VAP(J)-QuiD(J+l) 
DO  502  1=1, NC 

<■  h 


x  ( i ,  j+i  ):^;:)va(:,j) 

C  A  L  <_  S  u  3  P  T  (  T  E  M  P  (  J  *  1 )  ,  \  C  »  3  ?  E  R  R  ) 

G  J  I  '.  0  t-i  i  J  -  l  5  =  (j  .Q 

DO  5  3  4  I s 1 , N C 
504  OUlMCrit J^I)=QUlMOH( J»i)*7HAL?FCE\TnK( I ),ENTHL( I)#TrMP(J*l))*X(I# J* 
11  ) 

QC  =  VAP  C  J ) *  (  VAPKOH  ( J)  -QU  I  M3.1  (  J  +  1  i  ) 

GO  TO  443 
201  o=2 

LTIMESs: 

GO  TO  44 
22  i  3J.\3ALC  J)a(r  DThALP-VA?(  j)*VA?KOH{  J)-0UIM0H(J)*<QUID<  J)*RHDOT<J)  )*0 

i  a )  /  r  w  ( j ) 

DIFF =UER:VH(j)-QUMBAL(J) 
IF(A3SF( DIFF )-DELERR> 225,225, 222 

222  L  T I  X E S  =  L  T  I  M E S - 1 

IF (LT I K=S>224, 223, 223 

223  RrtO=RW(j) 
RWREF=RW( J) 

RW(J)=RW(J)+0.0i*DELTAU 
R/.»DOT(  J)  =0-01 

VAP( J)=QUID( J+l)-QUID(J)-RWDOT< J) 

DIFFO=DlFF 

GO  TO  44 

224  SLOPE =(DlFF-DlFFO)/(RW(J)-RWO) 
RWNEW=RW< J )-D IFF/SLOPE 
RWDOTC J)=(RWNEW-RWREF)/DELTAU 
VAP(J)=OUID( J*l)-OUIO( J)-RWDOT( J) 
DIFFO=DIFF 

Ra'0  =  RW(  j> 

RW  <  J>=RWNEW 
IF(LTIMES*3o)226,44,44 

226  PRINT  227.0U.MBAL(J)#DERIVK(J) 

227  F0RMAT(i3HREB0lLER  LOOP,2F20.6) 
GO  TG  30 

225  DO  226  1=1, NC 

X( i ,J-il=TX( I ,J) 

226  X( I ,J)=7X< 1/ J) 
253    NTALLY-,viTALLY  +  l 

IF(NTA^LY-NPRINT)4,17#17 

17  PRINT  1G#N'I  ALLY,DELTAU 

18  FORtfATCi'/HliNCREMENT  NO.  =  I5/34H0  INDEPENDENT  VARIABLE  INCREMENT  = 
1  FlQ.4i 

PRINT  19 

19  FORXAT(ll6H0ST6    KOL  FRA  SUM    H  UNBALANCE   H  DERIVATIVE    TEMPER 
1ATURE    LIQUID  FLOW   LIQ  ENTHALPY    VAPOR  FLOW   VAP  ENTHALPY       ) 

JSTG=0 
JLIM*NS*1 
DO  27  J=2,JLiM 

PRINT  2,-j,  JS7G,SUM(  J),CUNBAL(  J)  . DER I VH ( J > , TEMP < J ) , QU I D < J > , QU I MOH ( J ) 
1,  VAPC J)  , VAPMOH( J) 
2G  FORMAT ( J4,6Ei4.B) 
27  „STG  =  JS7G-rl 

PRINT  271,00,  QR 
27:  FCSXATdoKjCONDENSER  DUTY  =  Ei5 . 6/-7H0RE8O I LER  DUTY  =  El5.8) 

3G29  u=2»«ul« 

PRINT  25, JiTG, (X( I, J), I=1,NC),RW( J),RWDOT(J) 

25  FCRmaT(i9HC       STAGE  NO.  =  13/ ( 10F13 .10 ) ) 
29  jSI"G  =  JS7G-1 

6  7 


n?hint=npr:  v  r- 

IF(\T*lLY-NC\" 
30  STOP  77777 
END 


)  ^,3q,30 


SUBROUTINE  FOR  BUBBLE  POINTS 
SUBROUTINE  3U3?T(T,:_,3PERR) 

dimension  Gj;:x(:G;,A(io),S(io 
common  o  :dx, a,b,c 

EQUlLKF(A  -S,C,  V)=EXPF(A/(T-<6q 
K  T  I  M  E  S  =  i 

3  SUMY =Q.2 
DO  4  I=i/L 

Y  =  EQUlLKF(A(  I  ),3(  I  )  ,C<  I  ) . 7)*Qu 

4  SUMY=SUMYtY 
IF<A3SF(SUMY-1.C)-8PERR)6,6,5 

5  KTIMES=KTIMES-1 
IF(KTIKES)  7,6,6 

6  SUiviYO  =  SUrlY 
TO  =  T 

T  =  1+10.0 
GO  TO  o 

7  SLO?E=(SU,^Y-SUMYO)/(7-T  '•) 
TN=( (1.0-SUMY)/SLOPE)+T 
SUMYO=SUMY 

TO  =  T 


CO) 
*3+C*< T*460 .0 ) ) 


GO  TO  3 
6  tf  E  i  U  R  N 
END 
END 
OC  AND  OR  INITIAL  VALUES       11345.03732467 
GR  AT  CONVERGENCE       13247  -,23946193 
AlPH  BET     1.215010      .330>^C 
HOLDUP  RATIOS 

.000000  10.000000    .985907    .99*619 


13667.02693743 


999Q95   1.000000    .844210 


.86 


DROSRAM  DYNAMIC 
THIS  hROGRaM  HAS  BOTH  CONDFNSER  AND  PEBOILER  DELAYS  AND  ACCU 

DIMENSION  X(10»e50^.TX(in,5o).5<ir(10.50).Yir(in.5o),XISSF( 
ISSFdO.SO  ),ENTHK(10)  ,FNCM0).DlYV(50).DLYHV(50)  ,  DIYHL.C50) 
?  CMTHI  (1 0 )  .ENTHVM 0> »FNTHW(1 0) ,FK(1 0) . ALPHA(3) ,BETA(3) , 
3DFRIVM(lO),TEMP(5n) .VAP(5m,QU!D(50).FL(50).FV(50>.SSFl(5 
4n),SSPL(5n),SSPV(50)  .SUM(5fl).l.'A3MnH(50).OUIMOH(50).DERIVH 
5 A'  (50)  »HFL(50)  , HFV(50) .HSSFL  (50),HSSFV(50),RW(50),RWDOT<5 
6RHOQO(ln).CON(ln)'.FFen(lO).OUTOX(1S).A(lo)»B(lO>»C(lO)tDL 
7.nELAYL(5n).DELAYlH(50).DFLAYX(10.50).P(20»40) 

COMMON  OUIDX,A,B,C 

EOUILKF(A,3,C.T)=EXPF(A/(T*«60.C)*3«-C*<Tt-460.0)  ) 

THALPF(Y,2.T)=Y*T+Z 

RFAD  2,wC»NS, JFEED.NCNTROl 
2  F0RMATUI5) 

RFAD  333, (ALPHA(M,BETA(I).GAMMA(I),I=1,3) 
333  FORMATUF12.6) 
1  FORMAT(flFlO -4) 

JI.IMsNS+2 

READ  1,  (  (X(  I  , J)  .  I=l.NO  .  J  =  l, JLIM) 

RFAD  l»((XIF(I.J).I=l.NC).Jr2.NS) 

RFAD  1,((YIFM.J).I=1.NC).J  =  2.NS) 

RFAD  1. ( (XISSF( I, J) . 1=1. NO . J=2»NS) 

RFAD  1, ( (YISSF( I , J)  . I=l.NO . J  =  2.NS) 

RFAD  1.  (TEMP( J)  , J  =  l.  JL  IM) 

RFAD  1. (VAP( J), Jsl. JLIM) 

RFAD  1,  ( OUID( J)  . J=l  .  JL  IM) 

RFAD  1, (FL( J),J=2.NS) 

RFAD  'it  (FV(  J)  .  J  =  2.NS) 

RFAD  1, (SSFL( J) . J=2.NS) 

RFAD  1. (SSFV( J) ,J=2.NS) 

RFAD  1, (SSPLC J)  , J=2.NS) 

RFAD  1. (SSPV( J) . J=2.NS) 

RFAD  1, (VAPMOH(J) . J=l. Jl iMi 

RFAD  1,  (QiMMOHC  J)  .  J  =  l.  Jl  IM) 

RFAD  1. (HFL( J) . J=?»NS) 

RFAD  1. (HFV( J) , J=?,NS> 

RFAD  1.  fHSSFl.  (J)  .  J=2»NS) 

RPAD  1, (HSSFV( J) . J=2.NS) 

RFAD  1. <RW( J) , J=2. Jl IM) 

RFAD1, (RWDOT( J) . J=2.NS) 

RFAD  1, (ENTHK( I ) , 1=1, NO 

RFAD  1, (ENTHL( I ), 1=1. NO 

RFAD  1, (ENTHV( I ) , T=1 ,NO 

RFAD  1, (ENTHW(  I  )  ,  1=1. NO 

RFAD  1,  (A(  I  ),  1=1,  NO 

RFAD  1,  (B(  I  )  ,  1=1, NO 

RFAD  1,  (C(  I  )  ,  1=1,  NO 

RFAD  1,DELTAU,BPERR.DFLFRR 

RFAD  1, (RH0QD( I ) . T=1,NC) 


MULATOR 
10.50)  ,YI 

GAMMA (  3)  , 

0)  ,SSFV(5 
(5o),OUNB 
0). 
YX(10,50) 


DO  560  I  =1  ,50    «< 

DFLAYLCl  )=0UID(3) 
DFLAYLH(L)=QUIMOH(3) 


(=«? 


fed      \/W  Fccw  To    ccA/0£#Se£- 


UNO    virfb*    Fcevs  F.e'oM    &£86tc&e 


Dl  YV(L)=VAP(\SJ 

o  i  v  h  v  ( ;.  >  =  \i  a  ?  ?>  n  h  (  n  s ) 

Hi  "HI.  (L)=1U!K0H{!yiR*l)« 

D  n   -5  61    r  =  1 ,  \  c 

i)Fl.AYX(  T  ,  I.  )=X<  I  ,3)  *  ~ 

561     Dl   VX(  I  ,1    )rX(  !  ,NS*1  )* 

560    CONTIMJP 

*'FIR  =  4.A2 

H  W  =  n  .  2  5 

AREA=2l.49 

rn=n.2762A 

C  T  =  0  .  4  2 

GR=3?.? 

RH0L=a .  n 

DO    1  CO     1  =  1  ,NC 
100    SriCL  =  RHnL    *    X (I . JFEFD,*RHOQD( I , 
CnNST=CT*SGRTFC2.*GR)*tfFIR 
DFNOM=CnNST*RHOI 

HT  =  HW    *     (;T)/nFMnf-n**.6667*(CUlD(  JFEED)  )**.6667 
W I  RSS  =  RHOt_*AREA*HT 
A  i  ?  H  =  A  R  F  A  *  H  W  /  W  L  R  S  S 

3PT  =  ArycA/^LrJSSvfrn/C0i\'ST)**.6667 
I.A9FI.     =RHAI.PH    BF7 
P3IMT    6.LAnEl.,AI  PH,RE7 
6    FnRMAT(A8,9Fl2.6) 

nn  9  j=^,ns 

R H  0  L    =0.0 
DO     13     1=1, NC 
13    RHOL  =  RHDi.     *    X(  I  , J)*RH0QD(  I  ) 
9    RW  f J)=AI  PH*RHOL    *    GFT*(RH0L*-.3333)*(QUID( J)**. 6667) 

p  r  I  m  t   innn#(RW(.j>.j  =  i.ji  I  m> 

1000     C"0RMAT(15H    HOLDUP    RATIOS     ./.12F10.6) 
aH?     =     0.1-DbLTAU 
NTALLYsh 
\!  P  R  I  M  T  =  n 
4    JPATH=1 

JVARY=JFEFD 

ji  n*'  =  JFFEn 

JHIGH=NS 

GO  TO  444 
445  J!  0/.'  =  3 

JHl3rl  =  JFFFn-3 

JVARY=JWlGh 

N!JM  =  -1 
4  4  4  OH  ?  6  M  =  J  I.  C"; .-.  .  J  ri  I  G  H 

J  =  J  V  A  R  Y 

JVARY;  JVARY-NUK 

LTlMESs-i 

44  nn  -  -.   ist.i'.c 

t  x  < : ,  j  -  - )  =  * :  '■ .  j  - 1 ) 

r  x  ( ; .  j  )  =  x  ( : , . > ) 
ii   Tx(r,j-i>=x(i,j*i) 

on  -.  2   i=i  .\c 

3ii  =  =  0.0 

=  rJ:<sEOJTL.''F(A(  I  ),PM  )  .C(I).TEMP<j)) 

TFRX  =  VAP(  .;  :  < F O .< - 0 1 '  T  D  (  J  ) 

GT\(  ;)  =  (■•'■'./(.;- 1  )  *FGH  II  KF(A(  I  )  jBv  I  )  *CU  )  .TEMP(J-l)  >*TX(  I*  J-l)  )+OU  ID 

1  ( .  j » i  ;  »•  ~  v  ■'  ' ,  j  -  j  ) 

FF  =  n(  I  )=r"Lf.J)*XIFf  I  .J)*FV(J)»YIF(  I,  J)+SSFL(J)*XISSF(I,  J)+SSFV(J)*Y 
1  TSS"(  I,  .15 


GUTPIJT  =  TX'(  I  ,  J)*(SSPI   (J)+SSPV(J)*EQf<> 

7o 


SAYs  f  CO?.' (  I  )-TX(  J  ,.j)»r  TESM*S*-SO" 

n  n  5  k  =  -, .  a 

T  X  <  !  .  J )  s  r  •/  <  :  ,  j  3  ♦  G  A  M  M  A  f  K  3  »  (  C  A  v  -  r, 
aiJE  =  ALPWA(K)*CAY*RE'S'AfiO*Qiii= 
5    CAYsf CON< 7 )-TX( I ,J>*r 7PRM+RWD3T 
12    73C  U  ,  J  )  =  Tx  (  I  .  J  3  «-C  A  Y/6  .•  G  -QUE/3  .  0 

31   ssiM(J)  =  n.n 
on  7   1=1 ,nc 

7  StJM<  J)sSJJM(  JWTX(  T  .  J) 

DO    B     1=1 ,  MC 

7XU  ,  J)  =TX<  I  ,  J)  /SDKC  J3 

8  QIIIOXC  1  )=TX(  !  .J) 

CALL    HUB»T(TEMP( jl.MC.BPERR) 
VAPMOH(.l)=fl  .f, 
OU!MOH(.l)sO  .0 

nn  15    i=i,mc 

EK(  I)sEOl)!l.Kr(A(  I  )»B(  T  >  .C(  I)  .TE 

DPRIVM( I ) = (CnN( I ) -TXf I , J3*(FK( I 
1     TX(  I  .  J)  +  (SSPL(.J)+SSPV(.|)*FK(  I  3 

DHIMOKC  .PsOUIMOHC  J3+TPAI  .  PF(ENTH 
15  VAPMOHf J)=VAPMOH( J)*ThAI  PF(ENTH 
1  ) 

TFRM 1=0.0 

TPPM2=0  .  0 

TPRM3  =  0  .  G 

T  F  R  M  4  =  0  .  C 

n  n  ?  i   :  =  ? ,  n  c 

TFRS; =TFPMl*f ENTHK( I )*TFMP( J)+E 
TFRM?  =  TF-.:i?  +  FK(J3*DFRTVM(T) 
TPKX3  =  TPPy3  +  EK(T)*TXfr..l)*(A(n 

?!     TFPMi=TFPM4»FNTHK( I )*TV( I . J) 

nPPIVH(J)=(TFRMl+(TFRM2*TPRM4)/ 

FHTHALPsVAPf  J-l)*\/APMrHf  J-l  )+OM 

1  (J)»HF"V(J)*SSFL(  J^*HSPFI   ( J3+SSF 

33  OllNBAL(.i)  =  (FDTHALP-OUTMnH(  J)*(S 
1  (<?S=»VJ  J3-VAP(  J)  )  )/RW( J) 
DTFFrOERIVhC J)-QUNBAL(J) 
IF(  ABSFf  niFF3-DFLFRR3  25.25.22 

2?  l.TIMES  =  l  TiMES-1 
R^Ol  =n .  0 
DO  lol  T  -1  ,  r*:C 
101  RHOLsRHfk+TXt  I  ,  J)*RHnr?D(  I  3 
pyPrRHQl »*.333o 
IP>'LTIMFS)  24.23.23 

23  R  W  f)  =  R  W  ( . )  3 

R  '.-•  R  P  F  =  R  W  (  J  5 
5/;  (J)  =  RW(J3  -  ph; 
Ra'DOTC  J)r.i 

".;,!  D(  J)  =  (  (R'.-'f  J)-AI  PH*RHr»L)/(BET 
?:;MF  DS  =  VAP(  J-l  >+Ottin(w'*1  3+Fl   (J) 
v&3(  J)aSiiviFDS-QUin(  J)-RWDnT(  J) 
35    ";  ;F-C  =  D  I  FF 
r:.-;    T'l    4'. 

24  SI  (V->E=(n:FF-niFF(n/(R!«iC  J)-RWO) 
=>>:>  =y.'  =  RH(  J)-DIFF/Sl.nPF 
R'/iB:3T{  J3  =  (K/.'NEW-RWRFF3/nEl  TAU 

D 1 1 !  P.  f  J  3  =  '  {  R  w  N  E  'w  -  A I  P  H  •  P  H  n  L  3  /  i  ri  c  T 

\>i.->(  j)  =sii«FOS-cuin(.;3-Rwonv(w»> 

rj  T  F  F  n  =  'J  I  F  F 

»•#  D  =  RW(.I3 

-■  /: !  J  )  =  ?.  A '-;  F '.-. 

IPf  l.TlMPS»3r,3     37,44,44 


(.J)  3  *r:rO(I  3-307 PUT3*DELTAU/RW< J) 

;= ) 

C  J)  >     *r-£ED(  I  ) -OUTPUT  )#DELTAU/RW(  J) 


MP( J) 3 

)*VAP(J)*QUID( J)-RWDOT< J))*FEED<  I  3- 

3  3 /RW( J) 

K(  I  3.ENTHL(  I  ).TEMP(^j)  )#TX(  I,  J) 

V(  I  ) ,ENTHW( I ) .TEMP(J) )*TX( I. J)*EK(  I 


NTHL( I  3  3*DERIVM< I  ) 

/(  (TEMP(  J)-*460.0)**2)-C(  I  )  ) 

TERM33 

IDC  J*D*0UIM0HC  J*1)*FL<  J3*HFl(  J3-FV 

V(  J)*HSSf"V(  J) 

SPL ( J3*RWD0T(J)*QUID( J)  >-VAPMOH< J>* 


•EXP)  3**1.5 

■FV( J)*SSFL< J)-SSFV(J) 


-EXP3  )**1.5 


7/ 


37  :.  ABPL  =  6HPVD  TEST 

GOTO  3;LAbFL'J^"NPAi{J>-n^rV-(J>'0UID^^VAP<J).RW(j)#RWD0T(J) 

25  DO  ?6  1=1 ,HC 

26  X(I,J)=TX(I,J) 
JPATH=JPATh+l 

GO    TO(l99,19Q,2ol > . JPATW 

199  IF(MTALl  Y-5M220. 200,200 

200  Ja     NS 
L=5f) 

CnMDMQH=DLYHL(L) 
OnsDLYV(L)*(DLYHV(L)-DLYHi  (L) ) 

cn\'Lio  =  nij]D(j+i)*n,iiD(j,2) 

LTlNlESsi 

RWREF=RW( J-1 ) 

RwDOK  J  +  !  )  =  OLYV(L)-CON'LIQ 

RW(J+l)=RWREF+RWDOT(J+l)*nELTAj 

205  00    206     1=1 , NC 

206  TX(I,  J+i'  )=X(  I  ,  J  +  l) 
DO    ?07     1=1, NC 

QUF  =  C1  . 0 

CAY=(niYV(L)*DLYX(I  .  L > - < CONL I Q  +  RWDOT < J  +  l > > *TX < I , J*l > > *DELTAU/RW ( J+ 

nn    20rt    K=1 , 3 

TX( I, J  +  1  )=TX( I, J  +  1  )+GAMMA(K)+(CAY-QUE) 

OHE  =  ALPHA(K)*CAY  +  RETAfK)*niJE 

208  CAY=(DLYV(L)*DLYXn  .  L  )  -  (  CONL  I  Q+RWDOT  (  J  +  l  )  )  *TX  (  I  .  J  +  l  )  )  *DEl_T  AU/RW  (  J* 

207  TX(I,J*1)sTX(I,j*l)+CAY/6.0    -    QUE/3. 0 
SUM( J*1)=0.0 

DO    209     T  =1  , MC 

209  SUM(J  +  1)=SUM(  J+1)+Ty(  J  ,  j«.j  ) 

on  2in   i=i, nc 

TX(  I ,  J  +  i  )=TX(  I, J  +  l )/SUM(J  +  l) 

210  numx(  i  j=tx(  i .  j  +  i) 

CALL    RUP.PT(TFMP<J  +  l>,l\'C.BPERR> 
nu  I  MOH  ( .1*1  )  =  0  .0 

on   2ii    1=1 ,mc 
FKC(I)=PQUILKF(A(T).BM)«CM)»TEMP(J  +  1>) 

DPRIVM(T  )  =  (DLYV(L)*nLYXM.L)-(CONLIQ  +  RWDOT(J  +  l))*TX(I.J  +  l))/RW(J  +  l 

211  QUIMOHU  +  1  )=QUIM0H(J  +  1)+THALPF(ENTHK(I),ENTHLU),TEMP(J  +  1))*TX(I#J 

TPRY 1=0.0 

TPRv.2  =  o.C 

TPRM3=0.0 

TFR^  =  0.fj 

on  ?j.2  T  =  i  , NC 

TPRMisTFRMl  +(EMTHK( M*TEMP( J+1)+ENTHL( I ) )*DERIVM( I ) 

TPRV!2  =  TFRM2  +  EKC  (  I  )  *DFR  T  VM  (  I  ) 

TFR^3=TPPM3*EKC(I)*Ty(I.J+l)*(A(I)/((TEMP{J*l)+46  0.}**2>-C(I)) 

212  TPRM4=TFRM4*ENTHK(I)*TX(I.J+1) 
DPPIVH(J+1 )=(TERM1+(TFRM2*TERM4)/TERM3) 

QliN=JAL  (.J*1  )  =  (DLYVfL)*CONDMOH-(CONLIO  +  RWDOT(J  +  l)  )  *QU I MOH ( J  +  l )  )/RW(J 
1+1  ) 
216    DO    219     T  =  1  ,  N C 

219  X( I ,  J+1)=TX( ! ,  J  +  l) 

220  LVARY  =  50    -< - 

DO    500    1=1,49 


Ol  Y\/(LVARY)=DLYV(I  VARY-J  ) 
Hi  YHV<LVAWY)=DLYHVfl  VARY-1  ) 
DIYHL(LVARY)=DLYHI   (I.  VARY-1) 


7*- 


pec *y  /AJTcztf&L  (sc+z'i 
/a/   coA/peA/se*  irwe  t'tx. 


501 
500 


5n? 
503 


504 
201 


601 
60? 


V  A  R  Y  - 1  ) 


6n* 
603 


605 


60  6 


60  7 


6  0  6 


no  5o;.  r=i  ,kc 

Ol  VX  (  !  ,l  v  £RV  )  sDl.  Y) 

LVA^VslvaPV-i 

L  =  1 

Di  YV(L)=VAF<  J) 

01  YHV(L) =VA°MCH( J) 

m  YiL(L  )  =P.  .  f! 

DO    502     1=1, NC 

EK(  I)=EO!iIUr<F(A(  !)  .P(T)  .C(I),TEMP(J>) 

01  YXM  ,1    )=ErK(  I  )*X(  I.  J) 

DO    50  3     !=1,NC  \— ■ 

ouinxc  i  )=ni_Yx(  i  ,l) 

7PMPC0N=TFNP( J) 

CALL    BURPTdFMPCON.NC.BPERCn 

nn    5f)4     r  =1  ,  NC 

Di  YWl(L)=niYHl(l  )+THAI  PF  (  PNTHK  (  I  )  ,  ENTHL  <  I  )  ,  TEMPCON  )  *DL  YX  (  I  ,  I  )4 

GO  TO  445 

J  =  2 

L7IMES=1 

L  =  50 

00  602  1=1 ,NC 

TX( I , J)sX( I , J) 

DO  603  1=1, NC 

O!JE  =  0  •  0 

FOK=E«UlLKF(A( I ),R(I) ,C(I),TEH 

TPPM=VAP( J)*EOK+QUin(J) 

CON( I )=nELAYL(L)*PELAYX( I ,L) 

CAY=(COM(  I  )-TX(  I  ,.J)*(TERM  +  PWDO 

OO  604  K=1 ,3 

TX( I,J)=TX( I, J)*GAMNA(K)*(CAY- 

QUP=ALPHA(K)*CAY*RETA(K)*QUE 

CAv=(COM(I)-TX(!,.l)*fTFRM*PWOO 

TX<  I  ,  J)=TX<  I  ,  J)+CAY/6.'0  -  OUc/ 

SUM ( J)  =n  -  0 

D"  6  o  5  I  =  1  ,  N  C 

S'JM(  J>=R'JM.'  J)*TX(  1  .J) 

DO  60 0  1=1  ,  NC 

TX( i ,J)=TX( I . J)/SUM( J) 

oiii nx(  n  =TX(  ;»j) 

CALL  HURPT(TFMP( J) .NC.8PEPR) 
\'AHV|nH(  J)=0  .  0 
3UI  HOH<  .J  1  =C  .  0 
OO  607  1=1, NC 

EKM  )=EOlJILKF<A(  I),R(!).Cl!),T 
DPR  !Vh(  I  )  =  (CGN(  I)-TXU  ,.))#(  EK( 
(3ilIMCH(.J5  =QU!MOH(  J)f  THAI  PF(ENT 
>/APMOH(.I)=VAPMOH(  J)+TWAI  PFfENIT 
1  ) 

tpr*i=o.o 
TPR*i2  =  n.n 
TPR«3so.n 

TPRX4=C.O 

0  0    60  6     1=1  ,\'C 

TSRX1  =  TP-<S::  -  I  FNTHK(  !  )  *TPM?  I  J  )  + 

TPRy?  =  TP=>r-?-  =  K(  1  )*DFRTVfM  1) 

T  PR  V3  =  TFtfttS  -  if  K  {  I  )  *  T  X  U  , .» )  -  (  A  <  I 

<4  =  t PRr-: .; -f>;thk  (  l  )  *" x  (  I  .  j ) 
D  =  R  I  v  H  ( . ; )  =  ;  t  f  R  m  i  ♦  c  T  r  p.  y  2  -  T  R  R  y  4  ) 
FnTHALP  =  DFl.  AV[.H(L^*npi  AVL(i    ) 
"  :  I  \  -  A  L  ( .  I )  =  f  F  D  T  H  A  L  P  -  V'  A  P  (  .1 )  *  V  A  P  ••; 
1  P  )  /  R  w  i  J ) 

diff  =  :;£pivi-(j)-oijm6al(J) 


3(  J)  ) 

T( J)  ) )*DELTAU/RW(J) 

QUE) 

T(J) ) )*DELTAU/RW(J) 
3.0 


FMP( J)  ) 

I )*VAP( J)+OUID( J)^RWDOT(J) ) )/RW( J) 

HKfl), ENTHL (I), TEMP (J))*TX(I»J) 

HV( I ) ,ENTHW( I ) ,TEMP( J)  )*TX(I.J)*EK(I 


ENTHLf I ) )*DERIVM( I ) 

)/((TEMP(J)*460.0)**2)-C(I)) 

/TEPM3) 

OH( j)-QUIMOH( J)*(OUlD( J)*RWDOT( J) )*Q 


7* 


l  iajs .  Decay  iA/re&/*L 

*4     S"0     AZ>'l  . 


\r(  ARSF  (";;-•>"  '— D?:.FL'P-P25.?25,?22 

222  LT  !  MES  =  I  '  j>  r-c.-: 
IF(l T I XFS 5  224,223.223 

223  R#I0  =  RW(.J) 
R  W  R  -  F  =  R  W  (  J  ) 
RW( J)=RWf j)+o . 01»HEI  TAj 

=?wDOT(j)  =  n.oi 

VAPf  .J)  =nF|.  AYL  (L  )  -OUI  Df  J)  -Rir'DOT  (  J) 
0  J  p  F  0  =  [j  I  r  F 
3D    TO    6 n 1 

224  Si  OPE=fniFF-niFFO>/(RM.n-RWC) 
RWN  =  W  =  RW(J)-niFF/<5Ln?r 
R W DOT  (  J  )  =  (  P kM t w -R'.iRF F  ;  / DEI  T  A u 
VAP(J)=nELAYL<L)-OuTD(J>-RWDCTfJ) 
DTFFO=niFF 
RW0=RW(J) 
R  W  (  J  )  =  R  W  N  F  U 
IF(l  riM,FS  +  3P  )226,603  ,601 

226  PRINT    J??7,QUNBAL(  J)  .DFRTVHf  J) 

227  F0PMATM3HRFB0II.EP    IOHP.2F20.rt) 
GO    TO    3  0 

225  00    228    I  =1  ,  rsJC 
X(I,J-1)=TX(I»J) 

228  X(I,J)=TX(I»J) 

L  v  A  R  v  =  5  n     ^ 

DO    60  9    1=1,49 

DFLAYL  (I   VARY)=DFLAY|.  (I  VARv-1  ) 
DFLAYLHfLVARY)=DEl  AY[MI  VARY-1) 
on    61 o     T  =  1  ,NC 

61. 0     DF-l.  A  YX  <  T  ,  I.  V  AR  Y  )  =  DFl  A  Y X  (  T  ,  I   VARY-1) 
609    _VARY=LVARY-l 

L  =  l 

DFLAYLd   )=QUID<  J  +  j  ) 

DFLAVLH(L  )  =QUIMOH( J+t ) 

D  n    6 1  i     I  =  1  ,  N  C 

611.     DFI.AYX(  T  ,L)=X(  I  ,  J  +  l)«* 

253     \iTAli_Y  =  MTAl.l.Y  +  l 

IF(\TALI  V-\,:3RINT)4,17,17 

17  °P!MT     lfi#NTAI  LY.DFLTAl; 

18  FORMATd' 7H1  INCREMFMT    <vO. 
1     F 1  n  .  4  ) 

=>  R  !  \!  T     l  o 

19  "P.^tT  (  *,  1.6r-0F  1  G    MDl  FPA  SUM 
1ATURE    LI  QUIP  FLO 

JSTG=0 

DO  2  7  Jr2,JLiM 

=>R!  \T  20,jSTG,SUXfJ)-CUVflAL(J).DERIVH(J),TEMP(j),QulD(J)»QUIM0H(J) 
1  ,  v:3(  J)  .  v'APXnH(.J) 
2  0  F  0  R  M  A  T  (  T  <: ■ ,  n  ?  1  4  .  6  ) 

27  JST3=JST^ 

3  R  I  \'  T  2  7  j.  .  "  0  .  0  R 
271  FORMAT  (.'.  8,n,  ■C0.\|DF\'<5ER  DUTY  =  F  '  .,  .  6/1  7H0  REBO  I  LER  DUTY  =  6.15.8) 
J  S  T  0  "  0 

no  2 R  j  =  2  , ., ^  :  x 

a  }  t  v  -     2 *  .  J S  T G  »  (  X  (  •  »  J  )  .  I  =  1 .  MO  )  ,  R  W  (  J  )  ,  R W DO T  (  J  ) 

28  "ORMATM9H0       STAGF  NO.  =  I  3/  (  1  0  F  13  .  1  0  )  ) 
2  9  JSTG=JSTG-1 

3  R  Ml''  2?9.:LY\/(5u).ni_v-.'(5  0).DLYHL(50).DLYX(l,50),DLYX(2.50), 
1  "  .  v  y  (  3  ,  n  1  )  ,  n  I .  Y  X  (  4  .  5  n  )  .  T  F  M  o  O  O  N 
229  -I^A'idoHTOP  STAGE  WAR  =  .  F 1 2 . 6 , / . 1 6HT0P  STAGE  VPH  =  ,Fl2.6./.16H 
1T.03  STAGE  LOH  =  ,F12.6, /,17HT0=  STAGE  PROD  =  .  4Fl2  .  lo  ,  / ,  17HC0NDENS 

7^ 


I5/34H0IMDEPENDENT  VARIABLE  INCREMENT 


H  UNBALANCE   H  DERIVATIVE    TEMPER 
lO  FNTHAL3Y    VAPOR  FLOW   VAP  ENTHALPY       ) 


1ER  TFMP  =  .F12.R) 

moRINTsm?RINt*10 

jr<NTAU  Y-NCNTROL)  4,30.30 
30  STOP  7777 

END 
C 

C        SUBROUTINE  FOP  RURBI  E  POINTS 
C 

SUBROUTINE  RUBPTf T,L,BPFRR) 

DIMENSION  GUIDX(10).A(10),B(10>,C(10> 

COMMON  OUIDX.A,B,C 

POUTLKFf A.BtC.T)=FXPF(A/(T*46o .0)-B+C*<T*46o.O)> 

KTIMESsi 

3  SUMY=o.n 
DO  4  1=1 ,L 

Y  =  E0UILKF(A(I),B(I).CU).T)*QUIDX(I) 

4  SUMYrSUMY-rY 
ir(AB3F(SUf'Y-1.0)-BPERR)8.8,5 

5  KTIMES=KT!MES-1 
ir(KT IMPS)  7,6,6 

6  SUMYO=SUMY 
TOrT 

T=T»10.h 
GO  TO  3 

7  SI  OPF=(SUMY-SUMYO)/(T-TO) 
TN=( (l.n-SUMY)/SLOPF)+T 
S'iMYO  =  SUMY 

TO  =  T 
TsTN 
GO  TO  3 

8  RPTURN 
END 
END 

OC  AND  OR  INITIAL  VALUES       1  1R4:> .  03732467       1368  7.02693748 
OR  AT  CONVERGENCE       13247.23948193 
ALPH  BET     1.915010      .330550 
HOLDUP  RATIOS 

.OOOOnO  lO.nOOOno    .98-3907    .994619    .999095   1.000000    .844210    .864148 


k .  ;  •  ■  • 

APPgA/O/Xig 

COLUMN    PARAMETERS    FOR     HYDRAULIC    EQI\fS. 
FEED   RATE     =    250  gai/mjn 

COLUMN    DIAM.   =     6  FT 

ACTIVE    TRAY    AREA  »  76%     TOWER    AREA 

WErR    LENGTH    *    77%    OF    TOWER  DIAM. 

WEIR    HEIGHT  «  0.25  FT 

CALCULATED  VALUES 

/=    4.62  FT    WEIR    LENGTH 

A»  21.49    FT*  ACTIVE   AREA 

Vv-0-25  WEIR  HEIGHT 

F«  16.576  MOLE/Mlf/ 

C*   0.42  WEIR    EQN    CONSTANT 

'•■■..  •  1 

THESE  VALUES  ARE  THEN  USED  TO  CALCULATE  COLUMN  CONSTANT! 

9      LIQUID  HEI 
"J 


%     -Yw+     (f-Tpa.j)  /3(l.)   '3      LIQUID  HEIGHT 
7/  ^-/SgVft'         x  J'  ON  STAGE] 


W|r>ss  8      PcAfpi     HOLDUP    REFERENCE   STAGE 

(FEED  STAGE  USED) 


OC  3 


$ 


UWLR)88       !  . 

[•CONSTANTS    FOR  EQN  12 

A  F  J 


76 


APPENDIX  IV 

COMPUTER  FLOW  DIAGRAM  FOR  NUMERICAL  SOLUTION  OF  GENERAL 
STAGE  EQUATIONS 

Summary  of  equations  is  as  follows: 


Kij  =  yij/Xij  =  exPEAi/Tj  +  Bi  +  CiTjJ  (3) 


Pij  =  ?trl  (4) 

Xij   ■  R^   t(ZinPUt   •   *>Ut>Ut>mass    "   ViJ1  <7> 

#hj   "  R^T   ^^   ~   S°UtPUt>energy  "  K^  <8> 

"WJ   =  ^Lj  +  PPLj     Lj*  <9> 

0 

R...  =  (Einput  -  Eoutput)         ..  (10) 

wj               r  'mass  overall  v  ' 


hi  =  aiTj  +  bi  <12> 


h.  =  Eh.X  (15) 


j      n   ij 


hj  =  E(3iTj  l'   bi)Xij  +  ^AjU   V  U'),^    r    x  / 
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APPENDIX  V 

FORTRAN  60  PROGRAM  OF  LEAST  SQUARES  FIT  OF  MODEL  RESPONSE 
AND  A  SYNTHETIC  TRANSFER  FUNCTION 

First  program  is  for  a  two  parameter  fit  of  transfer  function 

.3s 


G(s) 


(Tis+1)(t2s+1) 


where  G(s)  is  the  transfer  function  of  C,  component  of  overhead 
product  for  a  step  in  C,  feed  composition. 

Second  program  is  for  a  three  parameter  fit  of  same  transfer 
function  with  an  arbitrary  gain  adjustment,  K,  i.e. 


Ke"-3S 
G(s)   = 


(TjS+IXtjjS+I) 
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PROGRAM  Lb  13  -  (jS*/G    0*4****7*1*   P-*t\ 


FOR 


THb  DIMENSIONED  ARRAYS  HAVE  THE  FOLLOWING  MEANING-- 

A(10)  CONiAlNS  THE  PARAMETERS  *'HICH  ARE  TO  BE  VARIEI 


C 

C   XC200).  XB(200>,  AND  XC:200)  CONTAIN  THE  I  NDEPENDbNT 'PARAME TERS 

c     pan2Set^s!AI^   1HE  n8S6RVEJ  VALUES  0F    'HE  ^ncu5nT<theAobserRved  depen 


C       INPJT 

r  ^IrlU^L'r"111    BE  RE^QDJCED  0,  O'JTPUI  USED  FOR  LABELING. 

C  SECOND  CARD  FORMA  ( 12.  I  3. I  2. I  3, E10 . 2 ) 

C  IKs  NUMBER  OF  PARAMETERS  TO  BE  VA-IED 

C  IS=  NUMBER  OF  POINTS 

C  NOF  s  FUNCTION  NUMBER  (SEE  ABOVE) 

C  NINP  =  NUMBER  OF  INDEPENDENT  PARAMETERS 

r  nrSJhc2Dic,nS!D  iS  A  CR1T6rI0N  FOR  CONVERGENCE.  IF  THE  RbLAIIVE  VALUE 

c     c0on;eHrge"eDisLreCa2hNeGdS  BY  LESS  ™AN  EDSIN   ,N  Tw0  ^ccessive1  I TERA^^NS. 

C       THIRD    CARD(S)     -     FORMAT     (5614.6)     -     CONTAIN     f0W    ESTIMATE     OF     THE    PARAMETERS. 

C      OBU^nTSum^  '    C°"Ul"    ™E    '~^    ^     '«    "10,    «R 

C       DEPESDEC^vARI^L0ErT(56l4•6,     '    ^  ^    U»'     °«    ™E    °BSERVED' 

C      ?!SJ¥    ;^ENDENTU;i5Ii^64"6)     *    C°NTA^    XlI)'    °R    ™E    — ^    OF    THE 
C 

c 
c 
c 
c 

DIMENSION  A<10>.  X(200).  XB(200).  XC(200),  Z(200).  FINCR(IO) 
1R(200).  E(10).U(10).RR(10),W(10)  HNLR(IO). 

1  READ    106 

READ    100.     IK,  Is.NOF.i^lNP.EPSlL 
IF( IR)    200.200 -10 
10    READ    101.     (A(l     .     I    =    i.ifl)  J 

READ    101.     iFINCH(l).     I     :    l.IR)        ^ **/>e*"4*»>r*L    Pa/urs    *Ci\    f£.%hMiU\ 

READ    102-     (U(  I  i,  1=1.  is)    «*t — — ^ 

read  102.    <kr(  i  >,  i=i,  is)  *« *»#jtj   ♦  *■    SfttPreo    LC+wer  P*i.y**m**Lt 

READ  102*  (ki  J  )  .  i:i.  IS)  4-  e»«tS7D/*r*c  t#*l*#rl 
PRINT  2»(U(I).I=l,IS),(RR(l).i:i,iS),(W(i).i=1,is) 

2  FORMA!  (16Hu   ..  R   ^   I NP JTS , / , < 7E12 . 6 . / . )    ) 

READ  103.  (x(i  ,  i  =  i,is)  + •*c*/«r$e*   i  •  /.  i . a/ 

DO  18  I«l. lb  * 

WW  =  X( I )  -  i.o 

Z( I )  =  0-0 

DO  19  J=  1.  IS 
18  CONTINUE(1>  *  MJ,*U(-)*"J"(>J>*-^  -«_  4n>/W**f»*    *#)"*«"* 

IF(NINP-2>  20.15.15 
15  READ  104.  iXB( I ),  I  =  1, IS  J 


R.  IS, NOP .NlNP.EPSIL) 


( A( I ) .   I   3  1. l«> 


PRINT 
PRINT 
PRINT 
PRINT 


i?(  i  >.  I  =  LIS) 


tX( t).  I  s  1. IS) 


LIS) 

<C.Z.F1NCR.EPSIL.NUITR.RRQ,N0F,R,E) 


IT(NINP-3)  20.36.16 
16  READ  105.  (XC( i ).  I  «  1. IS) 
20  PRINT  106 
PRINT  109 
PRINT  107 
PRINT  108 
PRINT  109 
PRINT  HO 
PRINT  109 
PRINT  Hi, 
PRINT  109 
PRINT  112 
PRINT  109 

PRINT  HI,  tFINCR(I).  I  =  1,  IR 
PRINT  109 
PRINT  lt3 
PRINT  109 
PRINT  H4, 

109 

115 

109 

116, 
PRINT  109 
JFCNINP-2)  30.25.25 

25  PRINT  117,  CXB  I  ),  I  =  LIS) 
PRINT  109 
IF(NINP-3)  30.^6,26 

26  PRINT  118.  IXC* I ),  I 
PRINT  109 

30  CALL  LEASI ( IP, iS.A.X.XB, 
PRINT  109 
PRINT  119 
PRINT  109 

PRINT  120,  <A<  i ),  I  =  1.  IR) 
PRINT  109 
PRINT  124 
PRINT  109 

PRINT  120.  <E< I ).  I  *    X, IR) 
PRINT  109 

PRINT  12L  NOI  iR,  RRJ 
PRINT  109 
PRINT  122 
PRINT  109 

PRINT  123.  IRC  i  ).  I  :    LIS) 
GO  TO  1 

100  FORMAT  ( 12. 13. 12.13.610.2) 

101  FORMAT  (5E14.6/ 

102  FORMAT  (5E14.6 

103  FORMA1  <5E14.6> 

104  FORMAT  (5E14.6* 

105  FORMAT  (5E14.6  ■ 

106  FORMAT  (BOH 
1 

107  FORMAT  (99H  NUMBER  Of    PARAMETERS      NJMBEP  UF  P01N1S      FUNCTION 
1  NUMBER     NL'MEER  OF  INDEP  PAMAM      E3SILON> 

108  FORMAT  (8X. i2.?4X. I3.l7x. I2,19> , I2.12A.bl0.2) 
10V  FORMAT  (1H  ) 

110  FORMA.  <26h  THfc  PARAMETERS  FED  IN  ARE* 
HI  FORMAi  (1H  .5E14.6) 

112  FORMAT  (55H  THE  INCREMENTS  FOR  OBTAINING  NUMERICAL  DERIVATIVES  ARE 
1) 

113  FORMAT  (34H  THE  OBSERVED  VALUES  TO  BE  FIT  ARE) 


114  FORMAT  (5E14.6) 

115  FORMAT  (31H  THfc  INDEPENDENT  PARAMETERS  ARE) 

116  FORMAT  (1H  ,5E14.6> 

117  FORMAT  (1H  .5E14.6) 

118  FORMA1  (1H  -5E14.6) 

119  FORMAT  (36H  THfc  BEST  VALJES  OF  THE  PARAMETERS  ARE) 

120  FORMAT  (1H  .7E17.1Q) 

121  FORMAT  (24H  NUhBER  OF  ITERATIONS  =  ,  I2.50H 
1ES  OF  THE  ERRORS  IS  NOW  =  ,E12.6> 

122  FORMAT  (53H  OBSERVED  VlNJS  CALCULATED  VALUES 

123  FORMAT  (1H  .9E12.5) 


HE  SUM  OF  THE  SQUAR 
THE  BEST  FI [  ARE) 


124  FORMAT  (66H  ES  IMATES  OF  THE  ERROR  IN  EACH  PARAMETER  ARE  (STANDARD 


200  STOP 
END 


SUBROUTINE  LEA^T (  I  R,  I  s . A . X , XB, XC . Z . F I  NCR , EPS  I L , NO  I TR . RRQ , NOF , R. E ) 
IR  =  NO.  OF  PARAMETERS,  IS  *  NO.  OF  POINTS,  A  IS  ARRAY  OF  PARAMETERS. 
X  IS  INDEPENDENT  VARIABlE(  I    DEPENDENT.   FINCR  IS  ARRAY  OF  INCREMENTS  OF 
PARAMETERS.   EPSIL  IS   -FRACTIONAL-   ERROR  CRITERION.   NOUR  IS  NO   OF 
ITERATIONS  REQUIRED  (UP  TO  10).   R^O  =  SUM  OF  SQUARES  OF  RESIDUALS) 

NOF  IS  NUMBER  OF   HE  FUNCTION  USED  IN  -EON". 
E  IS  THE  ARRAY  UF  fcSTIMATEQ  ERRORS  IN  THE  PARAMETERS 

DIMENSION  AilO.'.  X(2Q0),  XB(200).  XC(200).  2(200).  FINCR(IO).  R(20 
10),  0(200.10),  DT(10.200),  DEL'10),  DS(10).  DHK10,10).  DP(IO.IO), 

NOITR  s  0 
DO  20  I  =  1,  IS 
20  R(I)  =  U\,     -  fcQN(A,X(  I  ),XB(  I  )  .XC(  I  ),NOF) 

PRINT  4.  (R(I).  1=1, IS) 
25  IF  (N0ITR-9>  130,130.101 
130  NOITR  s  IM  0 1  I  R  «  1 
DO  220  J  =  1,  TR 
A(J)  =  A(J)*FINCR(J) 
DO  15  I  =  1,  IS 
15  D(I.J)  s  EUMA.X(  I  ),XB<  I  ),XC(  I  )  ,NOF) 
A(J)  =  A(j;-FINCR(J) 
220  CONTINUE 

DO  30  I  =  1.  IS 

CONST  s  EQM  A.  x(  I  )  ,X8(  I  ).XC(  I  )  .NOF) 
DO  30  J  =  1. IR 
30  D(I.J)  s  (D( I ,J)-C0NST)/FINCR( J) 
PRINT  3. ( <D( I.. i). jxi. IR) , 1  =  1, IS) 
3  F0RMAT(29HDERP,ATIVES  OF  EACH  P  ARAME  TER  ,  /  ,  (  4F1  0  .  5  ,  /  ,  )  ) 
DO  35  I  =  l.  IS 
1.  IR 
Di I, J) 
1,  IR 
1.  IR 
0-0 
1,  IS 

DP(  I  . J)*DT( I.K)*D(K, J) 
CALL  GAUSS3  (IR,1.00E-30,DP.DJI.KER) 
IF(KEH-i)  120,37,120 
37  DO  40  I  =  1.  IS 


PI (  J.K)*DHK.  i  ) 


DO  35  J 

35  DT( J. I) 
DO  36  I 
DO  36  J 
DP( I.J) 
DO  36  K 

36  DP( I, J) 


DO  38  J  =  1.  IR 

DS(J)  «  0-0 

DO  38  K  =  1,  IR 

38 

DSU)  =  DS(J)> 

DO  39  L  =  1.  IR 

3  9 

DT(L.  I  )  =  Db(L 

*0 

CONTINUE 

DO  110  I  -  1.  if* 
DEL ( I )  =  0.0 
DO  110  J  =  1.  »S 
110  DEL(I)  =  DFL(  I  ■•♦DTC  I  .-)••*<  J) 

DO  10  1  -  I.Ik 
10  A(  I  >  =  AM  )*DEL(  I  ) 

DO320  I  =  IMS 
320  RM>  =  ZM>  -  bON(A.<M).XB(  I  ) -XC(  D.NOf  ) 

PRINT  4,  <RM  >  ,  Irl,  IS) 
4  FORMAT  (16HDIFI-  \H    OdS-CALC  »  /  .  7F 10  .  5  ) 

PRO  s  0.0 

DO  50  I  =  1. IS 
50  RRJ  =  RR**Ri  I  ).«2 

CRES  =  A8SF (RRO-RRP)  -  E-JSIL*R^J> 

RRP  =  RRQ 

ir(CRES)  100.100.25 
101  FORHAf  (20H  CONVERGENCE  FAILURE) 

PRINT  101 

GO  TO  100 
120  PRINT  1001 
1001  FORMAT  (16H  MA  WlX  SINGULAR) 
100  FIS1  =  13-1 

DO  150  I  *  1,1* 

SUM  s  0.0 

DO    140    J    =    1. I> 
140    SUM    =    SUM*Di J. : )*«2 
150    EM)    =    SURIF<H'-Q/<SU1«FIU)> 

RETURN 

END 

SUBROUTINE  GAUSS3( ^, EP . A, X . KER >                                      000)0 
DIMENSION  AaO-10).  <(10.10) 

DO  1  1*1. N  00030 

DO  1  J=1«N  00040 

1  X{ I»J)*0.0  00050 
DO  2  K=1,N  00060 

2  X(K»K)si.o  00070 

10  DO  34  L'l.N  00080 
KP=0  00090 
Z«0.0  00100 
DO  12  KsL.K  00110 
IF(Z-ABSF(A(K.i  )))11.12,12  00120 

11  Z=A8SF(A(K,L) )  00130 
KP=K  00140 

12  CONTINUE  00150 
IF<L-KP)13-20.20  00160 

13  DO  14  J=L,N  00170 
ZrA(L.J)  00180 
A(L»  J)=A<Kf-  .  J)  00190 

14  A(K»»J)sZ  00200 
DO  15  Jal.N  00210 
ZsX(L.J)  00220 
X(L#  J)sX(Kf-.  J)  00230 

15  X(KP»J)=/  00240 
20  IF(ABSF(A(L.L) > -EP > 50 . 50 . 30                                            00250 

30  IF(L-N)31.34.34  00260 

31  LP1=L*1  00270 
DO  36  K  =  LP1.N,  00280 
IF( A(K.L> >32.36,32  00290 

32  RAMO  =  A(K.L  )/A»L.D  0030U 
DO  33  JsLPl.N  00310 

33  A(K, J)=A(K, J)-hATIO*A(L. J)  00320 


i5 


DO  35  Jsi,N 

35  X(K»J)sX(K,  J)-r-ATlO*X(L.  J) 

36  CONTINUE 
34  CONTINUE 

40  DO  43  1  =  1,  IV 
II=N-1-I 
DO  43  Jsi,N 
S  =  0.0 
IF< II-N)41.43.43 

41  IIP1=II*1 
DO  42  K=l IP1.N 

42  SsS*A( II»K)*X(K, J) 

43  X<  II. J)«(X(  I  I , J)-S)/A(  II.  I  I  ) 

RETURN  °°«*° 

H£IJKN  00470 

50  KER=2 
RETURN 
END 


00330 
00340 
00350 
00360 
00370 
00380 
00390 
00400 
00410 
00420 
00430 
00440 
00450 


00480 
00490 


FUNCTION  EQN(A- X.X8, <C.NJF) 
DIMENSION  A<10 

,X  =  xvini). 

|EQN=.0  52»EXHF(-3.«X)/((A(l)»X<-l.)*(A(2)«X  +  l.)*X«ll)6.)l 
RE'URN ~ — ; T X • 

END  fua/ctioa/   to  te  r*,r-T 

END 
U  RR       k       INPUTS 
.321390E-01  .189860E-01  .940000E-02  .463800E-02  .148000E-02  .530000E-04  . 0OOO0OE- 

.254460E-01  .129234E*00  .297Q77E+00  .500000E*00  .  702922E-00  .  870766E-00  .974554E- 

.647420E-01  .139853E^00  .190915E*00  .208980E*00  .190915E*00  .139853E+0Q  .647420E- 


Irf 


i>YN    UF    XFR    FCI    OF    n  *M    <»sXR(-3S>/< rS*l>(BS»l>     '0    i  UP    PROD    R6SP    '0    ^6D    *>E«T 

NUMBER    OF    PARAMEIfcRb  NjKtfE-«    OF    P'JlNTS  FUNCTION    NUMBER  VJM8ER    OF     INDE? 

2  7li 

iHfe    PARAMETERS    ••  bL»  ^1  k    ARE 

(O         (U) 

.J67000E*03  .410000E*02 
THE  INCREMENTS  H-*  OHTAJNUG  NJMERICAl  DERIVAJIVES  AkE 

•100000E*01     10^000E*01 

TME  "JBSEHVEU  VALltS   0  BE  -' I  I  ARE 

•  778981E-02    -161b92E-02    .591618E-03    .277002E-03    .148540E-03 
.8673645-04    -537G52E-04 

THE  INOEPENOENi  I-  AkAKETERS  ARE 

.lOOOOOE-Ol    -  ?  0 :  ■  0  0  0  E  -  J 1    .300000E-01    .40OQ00t*01    .5000006*01 
•  600000E-01     ■  70l'000fc*01 

DIFF     IN    08S-CALC 

.00013    .ooriui    .ooooo   -.00000   -.00000   -.  00000   -.00000 

DERIVAIIVES  OF  EACH  PARAMETg:< 

-.00002    -.00C05    -.00000    -.00002 

-.00000    -.00001    -.00000    -.00000 

-.00000    -.00000    -.00000    -.00000 

-.00000   -.00000 

OIF  F  IN  08S-CALC 

-.00001   -.00U01    .00000    .00000    -.00000   -.00000   -.00000 

DERIVATIVES  0F  EACH  PARAMEIch 

-.00002    -.00005    -.00000    -.00002 

-.00000    -.00001    -.00000    -.00000 

-.0  0000    -.00000    -.00000    -.00000 

-.00000    -.00000 
DIF  F  IN  08S-CALC 

.00000    -.00000 


.10000 
DEKIVAI  IVES 
-.00002 

-.00000     .00000 
OF  EALH  PARAMETER 
-.00005    -.00000 

.00000 
-.00002 

-.00000 

-.00000 

-.00001 

-.00000 

-.00000 

-.00000 

-.00000 

-.00000 

-.00000 

-  .  i.i  0  0  0  0 
OIF  F"  in  OBS- 
.00000 

-.00000 
■CALC 
-  .  0  0 1  U  0 

.00000 

.00000 

-.ooooo 

00000    -.00000 


<E   ^est   valued  0*    iHfc  pahakEiers  are|     two   p*d4*<eTZ+    R/f 

•  3441729345E*03   .  4584798486E-02    I        ^'  C    "* 


fc^TlMATE>  OF  THE  HWh.  B  IN  EACH  PARAM-E'ER  ARE  'STANDARD  DEVIATION) 


NUMBER  OF  IiERAIluNS  =   3     THE  SUM  OF  THE  SUJAWES  OF  THE  ERRORS  IS  NOW  =   .289( 

OBSERVED  MINUS  CALCULATED  VALUES  0^  i  HE  BES  i  FIT  A-iE 

.86381E-07  -.13422S.-05   .27661E-05   .14834E-05  -.77337E-06  -.24053E-05  -.33047E 
STOP 
TIME.   1  MINUTES  AND   12  SECONDS 


u 


program  Lb  -5  —  (J  P*****TW*  F-tr) 

C   MKlTg  FUNCTION  SobHROGHAM  T)  CALCULATE   hC  ruNCilUN  iO  Bfc  uEAb   SJUAREQ. 

C   A  IS  DIMfcNSIUNED  10  AMD  CONTAINS   *E  P  A-  A>4b  t  ERS  ,0  BE  VAniED  SuCM  AS  T0 

C   MINIMIZE  THE  SuM   F  ERRJRS  SQUARED  *.  *  B  .  AND  X'J  ARE  THE  i  H-itfc  INDEPENDENT 

C   PARAMETERS  tVtu  KAY  WISH  TJ  USE  ONLY  ¥  >  NO"   IS  i  Hfc  FUNCNuN  NUMBER  -  USED  FOR 

C   BRANCHING  TO  UiFK-UENT  *ARTS  OF  THE  EQM  >Jd3RUGRAM  WHEN  ^eVERA:.  FUNCTIONS 

C   ARE  TO  BE  Fll  (FO-  SEVERAL  JOBS  TO  BE  D^NE  »  . 

C   IHE  DMEiNSlONEU  A.-wAyS  HAVE  THE  FJLLOWlNG  MEANINU-- 

C      A(10)  CUNIAlNS  THE  PARAMETERS  ,>HiCH  AR£  TO  BE  VARIED 

C   X1200).  XB(20(i).  AND  XCC200)  CONTAIN  THE  INDEPENDENT  PARAMtlE"S. 

C   Z  (  200  >  CONTAIN  .HE  08SERVE0  VALUES  OF  'HE  FUNCTION  (THE  OBSE-vED  DEPrNDtNT 

C   PARAMETERS). 

C   riNCR'10)  CONTAIN   THE  MAGNITUDE  OF  THE  INCREMENTS  FOR  T Ht  ^AMAME'ERS  A(lO) 

C   SO  THAT  THE  PEUGRAM  MAY  TAKE  NUMERICAL  DERiVAIJvEb  WITH  hEAsONAHLE  ACCURACY. 

C   R<200>  CONTAINS  irtE  DIF"ERE^CE  BE • rfEEN  OBSERVED  U<200))  AND  CALCULATED 

C   fc(10)  CONTAINS  THE  ESTIMATES  OF  THE  ERRORS  Of  THE  PARAME'ERS  A<10>  AF  i  ER 

C   ITERATION. 

C       INPJT 

C   FIRST  CARD  -  ULL  BE  REPRODUCED  ON  OUT^i  uSEU  FOR  LABELING. 

C   SECOND  CARD  FORMA  (  12. I  3 . I  2 . I  3. E10 . 2 > 

C   IR  =  NUMBER  OF  PARAMETERS  To  BE  VAwIED 

C   IS  =  NUMBER  Ot-  ^OI.^TS 

C   NO*  «  FUNCTION  NUfHER  (SEE  ABUVE) 

C   NINP  =  NUMBER  i>   INDEPENDENT  PARAMETERS 

C   b^SIL  =  IS  USED  A-  A  CRITERION  FO-  CONVERGENCE.  IF  THE  RELATIVE  VALUE 

C   U"  THE  RESIDUAL  CHANGES  BY  lESS  THAN  EPS1N  IN  T*u  SUCCEbblVE  HERaTJOnS. 

C   CONVERGENCE  IS  REACHED. 

C    IHIRD  CARD(b)  -  FiWMAT  (5E14.6)  -  CONTAIN  ' OUR  ESTIMATE  Of      THE  PARAMETERS. 

C   AHO) 

C   FOURTH  CARD'S)  -  f- ORMA  I   (SEH.6)  -  CONiAlN  THE  INCREMENTS  FOR  I  HE  A(10>  ►  OR 

C   OBUINING  Nu^E^lCAL  DERIVATIVE. 

C   FIFTH  CARD(S)  -F 0^ MAT ( b  =  l4  .  6 )  -  CONTAIN  Hi),     OR  1  HE  OBSfcRVED. 

C   DEPENDENT  VARIABLES. 

C   SIXTH  CARD(S)  -  F'  *MAT  OE14  .6)  -  CONTAIN  X«I).  OR  THE  VAlUES  0^  THE 

C   FiRST  INDEPENDfcN;  VARIABLE. 

C   SEVENTH  CARDfb)  AID  FIGHT  CARD(S)  ARE  READ  IN  ONLY  If    MNP  IS  GREATER  THAN 

C  UNITY,  AND  CONiAlN  xB(I)  AND  XCU).  THE  SECuND  AND  THIRD  INDEPENDENT  VARIABLES 

C   STACKING  OF  JOBS  IS  PERMITTED  -  T*E  PROGRAM  STOPS  WHEN  Ji  READS  A  ZERO  FOR  JR 

C   (SECOND  CARD  INPli  ),  TH  JS  SEVERAL  BlANK  CARDS  AFTER  LAST  JOB  IS 

C   SUFFICIENT  !0  bTU-  ROUTINE. 

DIMENSION  A(10  ,  X(200).  XB<200>.  xC'200).  Z(200).  Fp.CR(lO). 
1R(200).  6  1 1  U  )  .  <  10  )  .RR'i  "  0)  .*(  10  > 
1  READ  106 

READ    100-     I".  I- .NOF. 'JIMP.EPSIL 

IF(  IR)    200-^00.10 
10    READ    101.     <  A(  l     .     i    =    1.  H> 

READ    101.     IMMUH(I).     I     =    1.  IR) 

READ    102.     «u(  1     ,  1  =  1. IS) 

READ    102.     (*&<:).  l*L  li> 

READ    102. (Wi I ) . |«1. IS) 

PRINT     2,(b(I)..=l.IS).(HR(I).I=l.lS).(WiI).I=l.IS) 
?    fQRMAI     (16ht        -R       W        IMP JTS./. < 7E12.6./. )        > 

t'EAD    103.     (  a(  1     .     I     =    LIS) 

DO    18     1  =  1.  lb 

-  W    =    X  (  I  )     -     1  .  c 

i  (  1  )     -    0-0 

00  19  J=  1.  lb 
U  H\)     -    Z  i  ■  *  -  (J>*0(-)«  <R<  J  )••*•* 
18  CONTINUE 

IF(NINP-2>  20.15,15 
1*5  «EAD  104.  i  XR(  )  .  I  ;  LIS) 

S7 


IFCNINP-3)  2  0.16.16 


I  =  1.  IS) 


1.  IR) 


1.  Ik 


'?<!).  I  =  LIS) 


X  (  I  )  ,   I 


1.  IS) 


16  READ  105-  uC(  i  ) 
20  PRINT  106 

PRINT  109 

PRINT  107 

PRINT  106,  (  IK.  IS.NOi-'.NIMP.EPSIL) 

PRINT  109 

PRINT  HO 

PRINT  109 

PRINT  HI, 

PRINT  109 

PRINT  112 

PRINT  109 

PRINT  Hi  ,  iFp-CRCDi  I 

PRINT  109 

PRINT  H3 

PRINT  109 

PRINT  114, 

PRINT  109 

PRINT  H5 

PRINT  109 

PRINT  lib, 

PRINT  109 

IF(NINP-2)  30. 25. 25 

25  PRINT  117,  IXB( I),  I  =  LIS) 
PRINT  109 

IF(NlNP-3)  30,26.26 

26  PRINT  US.  IXC<  I),  I  -     LIS) 
PRINT  109 

30  CALL  LEAS!  (  IP,  I  S  »  A  .  X  .  X8  ,  <C,  Z  ,  F  I  NCR.  EPS  1 1. .  NO  I  l  P  .  RRG  ,  NOT-  ,H,E5 
PRINT  109 
PRINT  H9 
PRINT  109 

PRINT  120.  ( A( l ).  I  =  LIP) 
PRINT  109 
PRINT  124 
PRINT  109 

PRINT  120,  ' E( i ) .  I  =  L IR) 
PRINT  109 

PRINT  121,  .nICjI  H,  RRJ 
PRINT  109 
DRINT  122 
PRINT  109 

PRINT  123,  <H< I ),  I  -  LIS) 
60  TO  1 

100  P0RMA1  (  12.  13.  i2.  I3.el0.2) 

101  FORMAi  (5E14.6 

102  FORMAT  (5E14.6* 

103  FOP^A'  (5E14.6. 

104  FORMAT  (5F14.6 

105  FORMAT  (5bl4.6 

106  FORMAT  (80H 

1  ) 

107  FORMA."  (99P  .NUMBER  OF  PARAMETERS 
1  NUMBER      NUMBER  OF  INDPP  PARAM 


NJMBEP    Of     POIN 
E=>MLON  ) 


FUNCTION 


11)8  fO«^Ai  (3*  .  i2.*4X,  13.17k.  I2.19*..12.12x.el0.2> 

1.  09  F  O  R  H  A  i  (  1 H  / 

110  FOp*AT  (26H  THr  PARAMETERS  FED  IN  A*E) 

111  FOHMA I  ( 1H  ,5b1 4 .6 ) 

112  FOPHA.  (55k  the  iNCRc^fcMrs  FOR  OBTAlNlNu  NUMERICAL  DERIVATIVES  ARE 
1  ) 

13  FORMA."  (34h  iHk  QBSERVfcD  VALUES  TO  BE  FIT  ARE) 


fi 


114  FO*^A'  (bfcll.ft 

115  ►■0«MAi  Cilh       t-e     INOfcJ£^iN       3AvAYrT 

116  FO^MAi  (lh     .5fi4.6) 

117  ro^MA.  (lh     .5b!4.6> 

118  f  0<<*Ai  (1H     .5bi4.6> 

119  rOK^Ai  {  3b»"     Thfc-     8f;S;      '/Ac   jgV     .IF     T.^ 

120  FOKMA>  (1.H     .7F17.10) 

121  FO^MAi  <?4H     \L'f"H6«     U-'     iii^AtJOO     J 
1ES    'Jf    THt    b»'»0-S    IS    *"rf    i     .612- 6) 

12*    FOKMA'  (53n    Ob    bBVED    v  1  \  JS    C»tCU-_A' 

123  FO^lAi  (1H    .9bl2.5) 

124  FOwMAi  (66H  fcS  IMArFi  Jt-  THfc'  E*<Oa 
1  DEVIATION)  > 

200    STO3 
ENO 


i2.:>0H  WE 

■0    Vi.JrS 


Hfc    Bfcsl 


;iCH     PA-AMEitW 


S  Ju4« 


ANJABL) 


SUB^OJTI  lb     LfcA     KH.  1>.A.X.X3.*C.?.«"  !  m  J' .fe^S  I  L  .  NOI  Trf  .  *-t(j,  \Dt-  .n.fc) 
I  *>     =     Mj.     OF     PAwAMtlEPS.      la     =     NC.       ■►     PDi'-j     S.     A     I-     A«R*V     J*     PA«<A«ETEBS. 
X     IS     IMDEPEMUbM        AR!ABlE<      '.     Dfc°E     U:M'.        ~  I  MO     IS     ABUa-     J*      INC    FMtNTS    OF 
PAWAMgl=Pb.       6  F-  >  J  I.     IS       -►UACMQNAl.-       E^n"<     -3 I f E* I JN .       NOJiH     IS    NO.     OF 
MtHATIONS    RtUfI«fcD    (U»    CC    10).       B-'Q    =    Sij.l    J^    SQ. AHES    0*~    ^ESIUualS) 


HE    F0  jc  r  I  )N   uSED    IN    -t3\- . 

STIHAfgU    ERRORS     IN     rnt    aa-UMEifc-vS 

.    X(2J0).    XB(200).    >;:<200>.    2(200).    rj*CH(10>.    -M20 

PT(ln.200).  DEl  10).  DS(IO).  Dfl<10.10>>  DP(10»lO)i 


NQf  ) 


NUt-      |S    N  )*Bfk    of 
b     IS    THr    aRRA.     \jf 
OIMENSIOi    1110 
10).    0(200.10 > . 
2E(10) 
N0I 1«    =    0 
DO    20    I         1 • IS 
20    P(I)     -    ZH)  t-uMC  A.  <(  i  )  .XB(  I  )  •  XC(  I  ) 

PRINT    4.     (r  .  I  )        1=1.  IS  > 
25    IF    (N)IT-.-V>    140.130.101 
130    NOITR    s    MO]  -\i     ■     J 
DO    220    J    =    L     !fi 
A(J)     =    Aijt'H'.C^J) 
'   ijj    DO    15    I     ■    LIS 
15    DU.Ji     s    ELMAX  (  !  )  .  <8U  >.«C(  L.NOF  : 
"A  (J)     =     lUi-fi'C^lJI 
22  0    CONTINUE 

U0    30     I     •     LIS 
CONST     =     t^v  A,      (  I  ) 
DO    30    J    -"    LIS 
3  0    D  <  I  .  J  >     -     (  Ln  I 
o^lNT     3,i(l i J 


4(  \   > ,XC(  I  )  .NJF  ) 


)-CJNSr ) /► I\CP( J) 
).J  =  L  I">.  1-1.  I  -> 


1  FQ8MA".  (29HLtft  I  .  ATH/ES  OF  EACH  -=  A  J  AM&  >  E  ■*  .  /  ,  (  4f-  10  .  5  .  /  .  )  > 


:  1 
-  1 
=  I 


DO  35  I 

DO  35  J 
So    0T( J.  I  ) 

DO  36  I 

DO  36  J 

DP( I .J) 

DO  36  K 
36  DP( I . J) 
CA.L  GA 

!F(KErt-l 
3/  00  40  I 

U0  38  J 

US(J)  = 

DO  38  K 
38  DS(J)  =  Dbi 

DO  39  L  =  1 
3V  DT'l.  I  )  -  i 
40  CONTINUE 


Is 


L  I* 

1-  In 

(•  •  0 

Lis 

by  (  1 

.3    ( 

120.37.120 

LIS 

L  Ik 

0 

1  .  ik 

i ./)  •►.  pi  (  j.<  »*u  r  (k.  i ) 


J)*l)T  (  I  .K  )  »D(  K.J) 
N.L0Je-3D.DP.U     I  ,Kt-<) 


Sv 


DO  110  I  =  1.  iM 

DEL(  I  )  =  0.0 

DO  110  J  -  1.   :S 
110  DEL( I )  =  UfeL( i  *DT ( 1 .„>•<< J) 

DO  10  I  =  1.  iw 
10  A(  I  )  -    Ail  >*DEi.  (  I  ) 

DO320  I  =  1» IS 
320  R(I>  =  Z(i)  -  bON(A.  <(  I  )  ,XB(  I  ).XC(  I  3  .' 

PRINT  4, <R( I ). 1=1, IS) 
4  FORMAT  (16HUIFF  IN  OrfS-C AlC . / . 7F 10 . 5 ) 

RRQ  =  0.0 

DO  50  I  -  1. lb 
5U  RRQ  =  RRu*M  ]  )*#2 

CRfcS  =  ABSMHRU'-RR?)  -  E^SIL'R^-' 

RRP  =  RRU 

IF(CRfcS)  100.100.25 
101  FORMAT  (20H  CU  VERGEMCE  FAILURE) 

PRINT  101 

GO  TO  100 
120  PRINT  1001 
1001  FORMAT  (16H  MATRIX  SINGULAR) 
100  FIS1  =  IS-1 

DO  150  I  =  1.  Ik 

SUM  =  0.0 

DO    140    J    -    1.  lb 
140    SUM    =    SUM-in  J,  i  )**2 
150    E(I)     =    SuHlHHhO/(SUM»FlJl)) 

RETURN 

END 


SUBROUTINE    GAUSS3<  M.fc^.A.X.KtR) 
DIMENSION  A* 10 -10) .  X(10.10) 
DO  1  1=1. N 
DO  1  Jsi.N 

1  X( I , J)=0-0 
DO  2  Ksi,\ 

2  X(K.K)*1.0 

10  DO  34  L=1.K 
KP=0 

2  =  0.0 

DO  12  K  =  L,N 

IF(Z-ABSF( AiK.i  )  )  )11.12. 12 

11  Z  =  ABSF( A( K.L  )  ) 
KP  =  K 

12  CONTINUE 
iF(L-KP)13. 20,^0 

13  DO  14  J=L.N 
7=A(l. J) 

A(L» J)=A(KF . J) 

14  A ( K P . J ) = I 
DO  15  Jsl,* 
ZsX(L.J) 

X (L, J; =X<  Kh  , J) 

15  X(K3,J):/ 

20  IF( ABoF( All .  L )  -EP)5J.30.30 
33  IF(L-N)31.34,3^ 
31  l_Pl  =  L*l 

DO  36  K  =  Lt-:  .  K 

I  F  <  A  <  K  ,  I  -  >  3  2  .  3  6  ,  3  2 
3  ?  «AUO  =  A(K,L'/A  L  .  L  ) 

uO  33  J  =  LH1  -N 

3  1   A(K,J>=A(K,j;--ATIO«A(L.J) 


00010 

00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 
00310 
00320 


Vc 


DO  35  Jsl.-.  lwjJl 

35  X(KiJ)sX(a.-)--AT|0«<(L.J)  0  0  J  4  L 

36  CONT INUE  G035C 
34  CONTINUE  C036C 

40  DO  43  Isl»*  U037L 
I I=M*1- 1  OOieC 
DO  43  J=1#N  003VC 
S=0.0  0Q4CC 
IF( | l-N)41.43.*3  U041U 

41  !IP1«1 !♦!  0042C 
DO  42  K=i  If-l.N  00430 

42  S=6*A(  UfKi*X(K»J)  0044C 

43  M  U»J)siXl  I  I,  )-S>/Ai  I  I .  I  1  )  00450 
KE*=1  H046C 
RE«J«N  0047C 

50  KEP»2  00480 
PEURN 

END  00490 
C 

HJNCUON    bL-<A.*,XB.XC.N)F> 
DIMENSION    A* 10  ■ 

X     -    X/IQO  . 

rEQN=A(3)«bXHF  j  -  3  .  •  X  )  /  ( t  A ( 1  )  • X ♦ 1  .  )»(A(2)*x*l.  )  «  »  •  1 0  0  .  >  f 

END 

U   HK   M   INPUTS 

.321390E-01  .189660b  01  .94000JE-02  46380Ub-02  .148U00E-02  .530000E-04  .000000b' 

•  254460E-01     .129234E    00     .29707/E-00  .500000fti0     . 702922E-0 0     .  870  7*6E*00     .974554b 

.647420E-01     -139853b    00     .1?0«15E»00  208980b*r<0     .190915b*00     .139853E-00     .647420b 


?/ 


bYN  OF  >FR  FC   Uh  F  .  PM  K«c^l-3>)/(  iS«l)i85*l) 


UP  PROD  HESP  iO  ^'EED  PERT 


NUMBER  OF  PARAMETERS 
3 


NJNdE^  OF  JOINTS 
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FUNC1 ION  NUM3EP 
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DUMBER  OF  iMDEt 
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1HE  PARAMETERS  FbU  |h  ARE 
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•IOOOOOE'01    -100000E+Q1    .100000E-02 
1HE  OBSERVED  VALUES   0  BE  r"II  ARE 

591618E-03 


.778981E-02 
.867364E-04 


161692E-02 
537052E-04 


!HE  INDEPENDEN!  PARAMETERS  ARE 


.100000E-01 
.600000E+01 


200000E*01 
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OH  F  IN  OBS-CALC 

.  i]  0  0  0  0    -.00000  .00000 
DEPJVAIIVES  OF  EACH  PARAMETER 
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-.00002     .03116  -.00000 

. 01132    -  .  00000  -  .  00000 

-.0  0000    -.00000  .00  28/ 

-.00000     .00171  -.00000 

.00110 
UIFF     IN    OBS-CALC 

-.10000    -.OOfOO  .00000 
DERIVATIVES  OF  PATH  PAPAMETeR 

-  .00002 


-.00002 
. 01100 

-  .00000 

-  .00000 
.00107 


.00005 
.030  23 

. ooooo 

.00000 
.00167 


.1433* 
-.00J00 
-.00000 

.00279 
-  .00000 


0\y  F  IN  OBS-CALC 

.00000    -  .  0  0  0 u  0 
DERIVATIVES  OF  EACH  PAWAMfc 
-.00002 


-  •  00002 
.01099 

-  .00000 


.0  000  5 
.03021 
.00000 
.00000 


00000 

E  1  ■=« 

.14529 

-.00000 

-.00000 
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300000E+01 

.00000 
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-  .00001 
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-.00000 
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-.00001 

.00515 
- .00000 
- .00000 
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-  .00001 
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400000E+01 


148540E-03 
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00000 


00001 
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- .00000 
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.00000 
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.00000 
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Ul*-  P  IN  OBS-CAlC 

.  rJ  0  0  0  0    -  .  0  o  r<  U  0 


-.00J0U 


OOJnj 


.00:00 


OO'JOO 


00000    -.00000    -.coooo 


1HE  BEST  VALUED  U» 
.3620037633E*03 


K-6  PAHA^fcTtRS  ARt    .      v 
<45852J789£*02   .  53619582 /9E -  0 1 


EiTlMAlES  Of  THE  bHR.  R  IN  cACh  PAHAMe'ER  A-t   STANOAhU  DEvUMON) 

.1248909550E-00   .  381598 76 74e-01   . 1455978 737E -  04 

NuMBhR  0^  Ii'ERAiIuNb  =   3     r  -IE  SJM  Of-'  THE  5-HAREj  Of  THE  kHRORS  IS  NOW  s   .282 

OBSERVED  11NUS  LALCui  ATEO  ^AlusS  Of   HE  Br-   ; I T  ARfc 

.310619-07  -.83954i--06   .24157E-05   .10739E'05  -.11070E-05  -.26591fc-05  -.34959 
STOP 
tIME.   1  MINUTEb  AND   15  icCUMOS 


7^ 


APPENDIX  VI 


DIMENSIONLESS  TIME,  t,  DEPENDENCE  ON  REAL  TIME,  9 


From  the  definition 


^a 


and  the  specific  column  parameters  given  in  Appendix  III,  the 
numerical  equivalence  of  equation  (1)  can  be  calculated.   The  calcu- 
lation is  as  follows: 

^  Lr'sS  ~  PLj  ^Tj  I  evaluated  at  4th  stage  from  '  ' 

initial  conditions 


Using  the  X.  .  from  Appendix  VII 


From  Appendix  III 


and 


(WLr)ss  =  4.4228  (moles) 


Therefore, 

t  =  3.748  9 


1  t  unit  -  0.267  (minutes) 
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(3) 


*ri  =  *j  {WT^J  ^U  -4144  <ft>  <4> 


APPENDIX  VII 
INITIAL  AND  FINAL  STEADY  STATE  VALUES  OF  COLUMN  VARIABLES 
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